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Abstract 

Women  have  a  1  in  8  chance  of  being  diagnosed  with  breast  cancer,  and  1  in  30  will 
die  of  this  disease  over  her  lifetime.  Although  mammography  is  the  current  procedure  of 
choice  for  breast  cancer  screening,  radiologists  fail  to  detect  cancer  in  10  to  30  percent 
of  patients  with  breast  cancer.  Also,  only  10  to  20  percent  of  patients  reffered  for  biopsy 
based  on  mammographic  findings  prove  to  have  cancer.  Further,  the  malignancies  missed 
by  the  radiologist  are  evident  in  two-thirds  of  the  mammograms  retrospectively.  To  increase 
sensitivity,  a  double  reading  has  been  suggested.  However,  the  explosion  in  the  number 
of  screening  mammograms  makes  this  option  unlikely.  Alternatively,  a  computer-aided 
diagnosis  (CADx)  system  may  act  as  a  “second  reader”  to  assist  the  radiologist  in  detecting 
and  diagnosing  lesions. 

Clustered  microcalcifications  are  one  of  the  earliest  indicators  of  breast  cancer,  and 
are  detected  only  by  mammography;  30  to  50  percent  of  nonpalpable  cancers  are  mam- 
mographically  visible  on  the  basis  of  microcalcifications  alone.  Furthermore,  for  early 
breast  cancers,  screening  studies  suggest  that  70  to  90  percent  were  detected  based  on 
micro  calcifications  alone.  This  research  proposes  the  following  methodology  for  clustered 
microcalcification  detection.  First,  preprocess  the  digitized  film  mammogram  to  reduce 
digitization  noise.  Second,  spatially  filter  the  image  with  a  difference  of  Gaussians  (DoG) 
kernel.  To  detect  potential  microcalcifications,  segment  the  filtered  image  using  global  and 
local  thresholding.  Next,  cluster  and  index  these  detections  into  regions  of  interest  (ROIs). 
Identify  ROIs  on  the  digitized  image  (or  hardcopy  printout)  for  final  radiologic  diagnosis. 
Finally,  to  improve  detection  rates,  globally  optimize  detection  parameters  using  a  genetic 
algorithm  (GA),  then  locally  optimize  using  the  simplex  method. 

The  database  of  56  digitized  (12  bit,  100  /um)  full-breast  (20x10  cm^)  film  mammo¬ 
grams  contained  63  biopsy-truthed  clustered  micro  calcification  ROIs  over  28  cases.  This 
technique  demonstrated  a  true  positive  (TP)  case  detection  rate  of  96.4  percent  (27/28), 
and  TP  ROI  (54/63)  and  TP  image  (48/56)  detection  rates  of  85.7  percent  with  5.75  false 
positives  (FPs)  per  fuU-breast  image. 
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Clustered  Microcalcification  Detection 
Using  Optimized  Difference  of  Gaussians 


I.  Introduction 

1.1  Background 

Mammography,  along  with  physical  examination,  is  the  current  procedure  of  choice 
for  breast  cancer  screening.  Screening  mammography  has  been  responsible  for  an  esti¬ 
mated  30  to  35  percent  reduction  in  breast  cancer  mortality  rates  (58).  However,  in  1996 
approximately  185,700  new  breast  cancer  cases  will  be  diagnosed  and  44,300  women  will 
die  from  this  disease  (47) .  Additionally,  women  have  a  1  in  8  chance  of  being  diagnosed, 
and  1  in  30  will  die  of  this  disease  over  her  lifetime  (60) . 

Although  radiographic  mammography  is  a  well-studied  and  standardized  methodol¬ 
ogy  (35),  for  10  to  30  percent  of  women  diagnosed  with  breast  cancer,  their  mammograms 
were  interpreted  as  negative  (21,  65).  Further,  the  malignancies  missed  by  radiologists  were 
evident  in  two-thirds  of  the  mammograms  retrospectively  (21,  65),  Missed  detections  may 
be  attributed  to  several  factors:  poor  image  quality,  improper  patient  positioning,  inac¬ 
curate  interpretation,  fibroglandular  tissue  obscuration  (57),  subtle  nature  of  radiographic 
findings,  eye  fatigue,  or  oversight  (21,  65). 

To  increase  sensitivity,  a  double  reading  has  been  suggested  (21,  65).  However, 
the  explosion  in  the  number  of  screening  mammograms  makes  this  option  unlikely  (55). 
Alternatively,  a  computer-aided  diagnosis  (CAD  or  CADx)  system  may  act  as  a  “second 
reader”  to  assist  the  radiologist  in  detecting  and  diagnosing  lesions  (21,  65). 

1.2  Computer-aided  diagnosis 

Researchers  at  the  Air  Force  Institute  of  Technology  (AFIT)  have  been  developing 
modifications  to  automatic  target  recognition  algorithms  used  in  airborne  reconnaissance 
sensors  (19,  51).  These  modifications  allow  the  technology  previously  used  to  identify 
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targets  in  sensor  imagery  to  be  used  to  identify  suspect  areas  in  mammography  films,  and 
to  render  a  recommendation  to  the  radiologist  as  to  the  pathology  (malignancy)  of  suspect 
areas.  Previous  AFIT  CADx  studies  have  resulted  in  detection  successes  ranging  from  77 
to  93  percent  (10,  33,  38,  48).  The  objective  of  the  AFIT  CADx  system  shown  in  Figure  1 
is  to  improve  radiological  detection  and  interpretation  of  lesions  in  mammograms.  Using 
features  extracted  from  detected  ROIs,  artificial-intelligence  techniques  can  then  be  used 
to  determine  the  possibility  of  malignancy  (21,  65). 


Figure  1.  CADx  System  Methodology:  After  film  digitization,  the  image  is  analyzed  by 
a  detection  algorithm.  Computer-detected  ROIs  axe  then  identified  for  the 
radiologists  to  compare  against  their  own  detections.  Next,  the  CADx  module 
provides  a  computer  diagnosis  for  the  radiologist  to  consider.  The  radiologist 
makes  the  final  diagnosis. 
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1.3  Problem  Statement 


The  objective  of  this  thesis  is  to  design  an  automated  microcalcification  detection 
system  to  be  used  as  an  aid  in  radiologic  mammogram  interpretation. 

1.4  Scope 

This  research  was  developed  to  aid  the  radiologist  in  mammogram  interpretation; 
the  radiologist  always  makes  the  final  diagnosis.  Also,  this  thesis  makes  no  attempt  to 
detect  or  classify  masses,  or  other  mammographic  densities. 

The  database  of  56  digitized  (12  bit,  100  pm)  full-breast  (20x10  cm^)  film  mam¬ 
mograms  contained  63  biopsy-truthed  clustered  microcalcification  ROIs  over  28  cases.  All 
mammograms  included  a  pathology  report  indicating  location  and  diagnosis  of  biopsied 
regions.  The  data  set  included  approximately  2  films  per  case  with  the  following  distri¬ 
bution:  1  case  with  4  films,  2  cases  with  only  1  film,  and  the  remaining  25  cases  with  2 
films. 

Full-breast  image  optimization  would  have  been  extremely  costly;  approximately  30 
to  60  minutes  per  digitized  film  on  a  Sun  Ultra  or  Sun  Sparc20.  Thus,  52  256x256  and 
9  512x512  subimages  were  extracted  from  these  images  to  facilitate  algorithmic  develop¬ 
ment.  The  relationship  between  full-breast  images,  subimages,  and  ROIs  is  illustrated  by 
example  in  Figure  2.  The  61  subimages  contained  all  63  biopsy-truthed  ROIs,  and  a  physi¬ 
cian  (27)  annotated  abnormalities  identified  in  corresponding  pathology  reports.  These 
subimages  considerably  reduced  algorithmic  development  times,  and  reduced  computer 
memory  requirements.  Additional  information  on  digitization  and  database  management 
is  located  in  Appendix  A. 
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Figure  2.  Digitized  (12  bit,  100  /lim)  full-breast  film  mammogram:  (a)  20x10  cm^ 
full-breast  image:  bounding  box  identifies  biopsy-truthed  region  of  interest 
(ROI);  (b)  256x256  pixel  subimage:  the  malignant  ROI  contained  a  cluster  of 
microcalcifications. 
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1 . 5  Methodology 

As  shown  in  Figure  3,  the  following  architecture  is  used  to  focus  radiologist  attention 
to  computer-detected  ROIs.  First,  the  digitized  film  mammogram  is  preprocessed  to  reduce 
digitization  noise.  Second,  the  image  is  spatially  filtered  with  an  optimized  difference  of 
Gaussians  (DoG)  kernel.  To  detect  potential  microcalcifications,  the  filtered  image  is 
segmented  using  global  and  local  thresholding.  These  detections  are  clustered  and  indexed 
into  ROIs.  Finally,  these  ROIs  are  identified  for  final  radiologic  diagnosis. 


Detection 
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Figure  3.  Overview  of  Cluster  Detection  Methodology:  t  controls  the  DoG  filter  used 
to  enhance  calcifications;  /,  fcjo,  khi,  and  N  control  the  detection  stage  of  the 
system;  fiCsmin  and  control  the  clustering  module. 
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L6  Overview 


Chapter  II  reviews  breast  cancer  and  related  CADx  research  efforts.  Chapter  III 
provides  an  overview  of  detection  methodology.  Chapter  IV  reviews  detection  results 
and  analyis.  Conclusions  are  drawn  in  Chapter  V.  Appendix  A  contains  digitization  and 
database  information  used  in  this  research.  Appendix  B  contains  Matlab  and  C  code  used 
in  this  research. 
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II.  Background 


This  chapter  reviews  breast  cancer  and  related  CADx  research.  Specifically,  this 
chapter  reviews  the  following  topics: 

•  Breast  cancer:  statistics,  mammographic  characterization,  radiologist  detection  and 

diagnosis  of  breast  lesions 

•  Related  CADx  research:  computer  detection  of  microcalcifications 

2.1  Breast  cancer 

Breast  cancer  results  from  failure  of  normal  regulation  of  cell  differentiation  and 
proliferation  (34).  According  to  Kopans  (34),  “Breast  cancer  kills  through  metastatic 
spread  and  compromise  of  the  function  of  other  organs.”  Since  many  cancers  require  time 
to  develop  the  ability  to  spread  successfully  to  other  organs  (34),  early  detection  is  the 
goal  in  screening  mammography  (60). 

2.1.1  Statistics.  Screening  mammography  has  been  responsible  for  an  estimated 
30  to  35  percent  reduction  in  breast  cancer  mortality  rates  (58).  However,  in  1996  approx¬ 
imately  185,700  new  breast  cancer  cases  will  be  diagnosed  and  44,300  women  will  die  from 
this  disease  (47) .  Additionally,  women  have  a  1  in  8  chance  of  being  diagnosed,  and  1  in 
30  will  die  of  this  disease  over  her  lifetime  (60).  Although  mammography  is  the  current 
procedure  of  choice  for  breast  cancer  screening,  radiologists  fail  to  detect  cancer  in  10  to 
30  percent  of  patients  with  breast  cancer  (21,  65).  Also,  only  10  to  20  percent  of  patients 
reffered  for  biopsy  based  on  mammographic  findings  prove  to  have  cancer  (55). 

2.1.2  Mammographic  characterization.  To  improve  radiologist  performance  and 
reduce  mortality  rates,  researchers  are  pursuing  computerized  analysis  to  aid  in  mammo¬ 
gram  interpretation  (55).  Radiologic  interpretation  of  mammograms  includes  detection 
and  diagnosis  of  suspicious  masses  and  clustered  microcalcifications.  However,  once  an  ab¬ 
normality  is  detected,  one  of  the  first  challenges  is  to  determine  whether  the  abnormality 
is  actually  present  (2,  16).  Numerous  overlapping  structures  in  the  breast  form  summation 
shadows  that  may  mimic  true  abnormalities  (8,  16). 
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2. 1.2.1  Masses.  A  mass  is  a  space-occupying  lesion  shown  in  two  different 
projections;  otherwise,  the  term  “density”  may  be  used  until  the  three  dimensional  features 
are  confirmed  (2,  16).  Densities  that  mimic  masses  include  pectoral  muscle,  nodular  breast 
tissue,  and  summation  shadows  (2).  Breast  masses  may  be  characterized  by  shape,  margins, 
density  relative  to  surrounding  tissue,  associated  findings,  and  location  (2).  According  to 
Adler  (2),  “Masses  are  three-dimensional,  distinct  from  the  surrounding  tissue,  and  most 
often  asymmetric  compared  with  the  contralateral  breast.”  Mammmographically  seen 
margins  are  one  of  the  major  determinants  of  benign  or  malignant  status  because  masses 
with  indistinct  margins  have  a  higher  risk  for  malignancy  (2,  16).  Additional  features  used 
for  interpretation  of  masses  are  discussed  in  more  detail  in  Adler’s  (2)  and  Dorsi’s  (16) 
articles. 


2. 1.2. 2  Microcalcifications.  Microcalcifications  are  calcium  deposits  typi¬ 
cally  on  the  order  of  0.1  to  0.3  mm  in  diameter  (8,  56)  and  are  detected  only  by  mam¬ 
mography  (20).  Recognition  of  their  presence  is  the  challenge  for  radiologists  because  30 
to  50  percent  of  nonpalpable  cancers  are  visible  on  the  basis  of  microcalcifications  alone 
(20).  Furthermore,  for  early  breast  cancers,  screening  studies  suggest  that  70  to  90  percent 
were  detected  based  on  microcalcifications  alone  (20).  (The  terms  micro  calcification  and 
calcification  will  now  be  used  interchangeably  to  simplify  discussion  of  characteristics.) 

As  mentioned  earlier,  one  of  the  first  steps  in  radiologic  analysis  of  calcifications  is 
to  determine  whether  they  have  been  detected.  Mimics  or  pseudocalcifications  include 
deodorant,  talcum  powder,  pickoff  of  film  emulsion,  dust  particles,  and  scratches  (20). 
After  false  detections  are  eliminated,  calcifications  are  analyzed  for  potentially  benign  or 
malignant  features. 

Feig  (20)  lists  the  following  examples  of  potentially  benign  calcifications: 

•  “Eggshell  or  rim  calcifications”  typically  have  a  round  calcific  rim  surrounding  a 
hollow  sphere. 

•  “Spherical  or  lucent-centered  calcifications”  have  thicker  walls  than  eggshell,  have 
smooth  surfaces,  and  are  round  or  oval. 
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•  “Vascular  calcifications”  typically  appear  in  parallel  producing  a  “railroad  track” 
appearance.  These  may  be  difficult  to  distinguish  from  potentially  malignant  fine 
linear  calcifications  if  one  wall  of  the  vessel  is  calcified,  or  if  the  segment  is  short  or 
discontinuous. 

•  “Large  rod-like  calcifications”  follow  the  course  of  the  ducts,  radiating  from  the 
retroareolar  area  and  are  frequently  bilateral. 

•  “Coarse  or  popcorn-like  calcifications”  appear  as  a  circumscribed  mass  with  coarse, 
mostly  round  calcifications. 

•  “Milk  of  calcium  calcifications”  are  layers  of  calcium  in  the  dependent  portion  of 
cysts,  and  may  be  meniscus  shaped  when  seen  on  the  mediolateral  (ML)  view  and 
less  distinct  on  the  craniocaudal  (CC)  view. 

•  “Skin  calcifications”  typically  appear  as  1  to  2  mm  lucent  centered  or  solid  spheres. 

•  “Dystrophic  calcifications”  have  a  bizarre  shape,  typically  result  from  surgery  or 
radiation  therapy,  and  may  be  difficult  to  distinguish  from  potentially  malignant 
calcifications. 

Shape,  distribution,  size,  and  contour  are  the  major  criteria  used  to  distinguish 
malignant  from  benign  calcifications  (20).  Feig  (20)  and  Shaw  de  Paredes  (56)  list  the 
following  characteristics  of  potentially  malignant  calcifications  (assuming  one  or  more  are 
present): 

•  “Linear  shaped  calcifications”  are  typically  considered  suspicious;  however,  some 
benign  conditions  may  also  produce  similar  calcifications  (20). 

•  “Linearly  distributed  or  segmentally  distributed  calcifications”  are  extremely  suspi¬ 
cious  since  most  malignant  calcifications  form  within  ducts,  and  over  90  percent  of 
invasive  cancers  arise  from  ducts  (20). 

•  “Markedly  clustered  and  unilateral  calcifications”  are  also  extremely  suspicious  since 
malignant  calcifications  usually  appear  with  this  feature  (20).  The  greater  the  num¬ 
ber  of  calcifications  within  an  area,  the  more  suspicious  (56).  As  Shaw  de  Paredes 
notes,  “malignant  calcifications  tend  to  occur  in  tight  clusters  of  1  cm  diameter  or 


9 


less”  (56).  However,  Shaw  de  Paredes  cautions:  “...the  disease  and  the  calcifications 
can  be  extensive,  involving  an  entire  quadrant  or  an  entire  breast”  (56).  Finally, 
solitary  clusters  are  typically  considered  for  biopsy  if  no  similar  calcifications  exist 
elsewhere  in  the  breast,  or  in  the  opposite  breast  (20). 

2.1.3  Radiologist  detection  and  diagnosis  of  breast  lesions.  Although  radiographic 
mammography  is  a  well-studied  and  standardized  methodology  (35),  for  10  to  30  percent 
of  women  diagnosed  with  breast  cancer,  their  mammograms  were  interpreted  as  negative 
(21,  65).  Further,  the  malignancies  missed  by  radiologists  were  evident  in  two-thirds 
of  the  mammograms,  retrospectively  (21,  65).  Missed  detections  may  be  attributed  to 
several  factors:  poor  image  quality,  improper  patient  positioning,  inaccurate  interpretation, 
fibroglandular  tissue  obscuration  (57),  subtle  nature  of  radiographic  findings,  eye  fatigue, 
or  oversight  (21,  65). 

To  increase  sensitivity,  a  double  reading  has  been  suggested  (21,  65).  However, 
the  explosion  in  the  number  of  screening  mammograms  makes  this  option  unlikely  (55). 
Alternatively,  a  computer-aided  diagnosis  (CADx)  system  may  act  as  a  “second  reader” 
to  assist  the  radiologist  in  detecting  lesions  and  making  diagnostic  decisions  (21,  65). 

2.2  Related  CADx  research:  clustered  microcalcification  detection 

To  improve  radiologist  performance,  researchers  are  pursuing  computerized  analysis 
to  aid  in  mammogram  interpretation  (55).  As  mentioned  earlier,  one  research  area  of 
interest  is  detection  and  diagnosis  of  microcalcifications.  Although  researchers  generally 
desire  the  same  goal,  several  methodologies  exist  for  each  stage  of  the  detection  and/ or 
diagnosis  process. 

2.2.1  Database.  Several  research  groups  are  currently  investigating  detection  of 
clustered  microcalcifications  (9,  44,  46,  64,  68).  Most  investigators  use  different  databases; 
thus,  it  is  not  possible  to  compare  results  directly  (44).  Additionally,  no  universally  ac¬ 
cepted  criteria  exists  for  the  number  of  bits  required  to  represent  available  gray  levels  in 
films  or  minimum  resolution  required  for  detection  of  clustered  microcalcifications.  How- 
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ever,  several  researchers  use  a  sample  spacing  of  100/xm  per  pixel,  with  a  gray-level  reso¬ 
lution  of  10  to  12  bits  per  pixel  (9,  18,  44,  64,  68). 

There  is  also  a  wide  variation  in  the  area  of  mammograms  used  for  analysis.  Image 
area  is  an  important  factor  in  increasing  the  false  positive  detection  rate  since  the  larger  the 
test  area,  “...the  greater  the  chance  for  producing  false  positives  in  general”  (18).  Typical 
mammogram  image  areas  generally  range  from  8x10  cm=^  to  approximately  10x10  cm^ 
(9,  18,  44,  68);  however,  robust  detection  algorithms  should  not  be  limited  by  image  area. 
As  mentioned  in  Chapter  I,  we  use  approximately  20x10  cm^  (2KxlK  pixel)  full-breast 
images. 

In  the  literature,  several  CADx  researchers  (9,  44,  46,  68,  64)  develop  and  test  de¬ 
tection  techniques  on  the  same  database  (not  to  be  confused  with  a  common  database). 
Nishikawa  et  al.  (45)  have  shown  “...the  choice  of  clinical  cases  used  to  train  and  test 
a  computer-aided  diagnosis  (CAD)  scheme  can  affect  the  test  results  (i.e.  error  rate).’ 
Nishikawa  et  al.  (45)  also  argue,  “Because  of  the  strong  dependence  of  measured  perfor¬ 
mance  on  the  testing  database,  it  is  difficult  to  estimate  reliably  the  accuracy  of  a  CAD 
scheme.  Furthermore,  it  is  questionable  to  compare  different  CAD  schemes  when  different 
cases  are  used  for  testing.” 

Although  a  common  database  has  been  suggested,  this  would  require  standardiza¬ 
tion  of  pixel  size  and  number  of  gray  levels  used  (45).  Also,  if  the  common  database 
were  enlarged,  all  previously  evaluated  CADx  techniques  would  need  to  be  reevaluated  if 
images  from  new  digitized  imaging  or  film-screen  systems  were  added  (45).  Furthermore, 
Nishikawa  et  al  argue,  “...because  there  is  no  quantitative  or  rigorous  means  of  determin¬ 
ing  whether  a  particular  common  database  is  representative  of  the  clinical  population,  it 
is  difficult  to  predict  the  true  clinical  impact  of  a  CAD  scheme.” 

2.2.2  Cluster  detection  criteria.  No  universally  accepted  criteria  exist  for  scoring 
detections  of  clustered  microcalcifications  (63,  64).  Nishikawa  et  al  (44)  rely  on  visual 
interpretation  for  determining  true  positive  (TP)  detections  since  “...it  is  not  possible  to 
determine  the  x—y  location  of  micro  calcifications  that  are  not  visible  in  the  mammogram.” 
Nishikawa  et  al  do  not  define  their  criteria  for  false  positive  (FP)  detections. 
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Strickland  et  al  (63,  64)  rely  on  the  cluster  detection  criteria  proposed  by  Karsse- 
meijer  (32).  “A  cluster  is  considered  detected  if  two  or  more  microcalcifications  are  found 
in  the  region  of  film  (truth  circle)  identified  by  the  radiologist”  (63,  64).  Strickland  et  al. 
continue,  “...a  false  positive  (FP)  is  counted  if  two  or  more  erroneous  detections  are  made 
within  an  empty,  closed  region  of  0.5  cm  in  width  (63,  64).  As  Strickland  et  al.  note,  “One 
weakness  of  this  metric  is  that  at  high  sensitivity  thresholds  the  number  of  false  positives 
may  actually  drop  due  to  merging  of  previously  separated  false  clusters”  (64).  However, 
this  usually  occurs  at  FP  rates  too  high  to  be  of  any  clinical  value  (64). 

Chan  et  al.  (9)  score  a  signal  detected  if  it  was  within  0.5  mm  of  a  true  microcalcifi¬ 
cation.  Further,  a  cluster  is  scored  as  a  TP  if  its  centroid  is  within  a  cluster  radius  of  5  mm 
from  the  centroid  of  a  true  cluster  and  at  least  two  of  its  member  microcalcifications  are 
scored  as  a  TP.  Once  a  microcalcification  or  cluster  is  considered  detected,  it  is  eliminated 
from  further  matching. 

2.2.S  Noise  reduction.  In  screen-film  mammography,  the  film  is  used  for  detec¬ 
tion,  storage,  and  display  (15).  Since  full-area  detectors  suitable  for  digital  mammography 
are  still  in  the  experimental  stage,  an  alternative  way  of  acquiring  mammograms  in  digital 
form  is  to  digitize  a  film  mammogram  (67).  Although  the  main  limitation  in  image  quality 
should  be  the  granularity  of  the  film  emulsion  (67),  noise  is  introduced  from  the  process 
of  digitization  (15).  This  noise  may  later  be  detected  as  a  pseudocalcification. 

Although  noise  reduction  may  be  performed  prior  to  processing,  Nishikawa  et  al.  (44) 
and  Yoshida  et  al.  (68)  apply  a  noise  reduction  technique  during  the  detection  process. 
To  reduce  the  effects  of  noise,  they  propose  use  of  morphological  erosion  to  eliminate 
detections  less  than  3  pixels  in  area,  since  very  small  signals  are  likely  to  be  due  to  random 
noise,  and  not  by  the  presence  of  true  calcifications. 

2.2.4  Calcification  enhancement.  To  enhance  calcifications,  several  CAD  micro¬ 
calcification  detection  methods  employ  a  difference  image  technique;  subtraction  of  a  signal 
suppressed  image  from  a  signal  enhanced  image  (9,  18,  44).  In  practice,  the  spatial  kernels 
used  for  enhancement  and  suppression  are  combined  into  a  single  linear  filter,  equivalent 
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to  bandpass  filtering  (18).  Several  of  these  methods  use  a  box-rim  based  kernel  to  enhance 
calcifications  (9,  18,  44). 

Strickland  et  al.  (63,  64)  and  Yoshida  et  al.  (68)  use  wavelets  to  enhance  calcifica¬ 
tions.  Strickland  et  al.  (63,  64)  use  sub-band  decomposition  to  obtain  individual  wavelet 
matched  filter  responses.  Yoshida  et  al.  (68)  use  wavelet  techniques  to  decompose,  mod¬ 
ify,  then  reconstruct  wavelet-processed  mammograms.  Finally,  other  techniques  which 
have  been  successfully  applied  use  locally  adaptive  or  region  growing  methods  to  enhance 
contrast  in  film  mammograms  (12,  13,  14,  42,  66). 

2.2.5  Detection  of  signals.  Once  an  image  has  been  preprocessed  to  enhance 
differences  between  targets  and  background,  the  image  is  more  easily  segmented.  Several 
researchers  threshold  their  filtered  images  to  detect  potential  calcifications  (44,  64,  68). 
However,  thresholding  schemes  typically  vary  according  to  the  methodology  used  for  en¬ 
hancement. 

Nishikawa  et  al.  (44)  and  Yoshida  et  al.  (68)  use  adaptive  thresholding  based  on 
local  statistics  of  gray-level  values  within  a  51x51  mm  area.  Pixels  are  kept  if  their  value 
exceeds  2  to  3.8  times  the  local  standard  deviation  plus  the  local  mean  (44,  68).  Nishikawa 
et  al.  (44)  use  a  64x64  mm  region  centered  on  each  signal  to  reduce  false  positives.  If 
the  first  moment  of  the  averaged  power  spectrum  exceeds  3  mm“^,  the  signal  is  discarded. 
Yoshida  et  al.  (68)  use  texture  analysis  to  reduce  false  detections. 

Since  Strickland  et  al.  (63,  64)  generate  multiple  images  during  processing,  they  pool 
individually  thresholded  images  to  establish  a  test  statistic  at  each  pixel  location.  Each 
image  is  thresholded  at  some  fixed  percentile  of  the  corresponding  histogram  (the  authors 
do  not  report  these  corresponding  percentiles).  These  images  are  then  summed  pixel- wise. 
Pixels  with  values  of  3  or  more  are  passed  to  the  final  detection  map. 

2.2.6  Detection  of  clusters.  Recall,  Feig  (20)  listed  “markedly  clustered  and  uni¬ 
lateral  calcifications”  as  extremely  suspicious  since  malignant  calcifications  usually  appear 
with  this  feature.  Chan  et  al.  (9)  and  Nishikawa  et  al  (44,  46)  also  argue  true  microcalcifi- 
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cations  of  clinical  interest  typically  appear  in  clusters.  Thus,  these  investigators  use  fairly 
similar  methodologies. 

Chan  et  al  (9)  use  a  1  cm  diameter  neighborhood  criteria;  if  the  number  of  signals 
found  within  the  neighborhood  is  greater  than  an  input  minimum  number,  the  signal  is 
kept.  Similarly,  Nishikawa  et  al  (44,  46)  cluster  detected  signals  by  passage  of  a  32x32 
mm  box  over  the  entire  image.  If  3  or  more  signals  exist  within  the  box,  the  box  is  passed 
to  the  output  image  as  a  ROL 

2.2.7  Calcification  detection  results.  Chan  et  al  (9)  evaluate  the  success  of 
their  convolution  neural  network  at  detecting  clustered  microcalcifications  over  52  100/^m, 
12-bit  films  (size/area  not  reported).  To  demonstrate  the  dependance  of  the  performance 
of  their  CAD  scheme  on  the  database,  they  (subjectively)  divide  their  films  into  obvious, 
average-subtle,  and  subtle  sets.  For  the  obvious  set,  they  report  a  TP  rate  of  100  percent 
with  less  than  0.1  FPs  per  image.  For  the  average-subtle  set,  they  achieve  a  93  percent 
TP  rate  with  one  FP  per  image.  Finally,  for  the  subtle  set,  they  report  a  TP  rate  of  87 
percent  with  1.5  FPs  per  image. 

Nishikawa  et  al  and  Yoshida  et  al  use  the  same  database,  as  Yoshida’s  article 
is  an  earlier  investigation  of  a  wavelet  enhancement  technique  (44,  68).  Yoshida  et  al 
achieve  a  detection  rate  of  85  percent  with  5  FPs  per  image  over  39  lOO/xm,  800x1000 
pixel  (8x10  cm^),  10-bit  images  (68).  Similarly,  over  this  identical  dataset,  Nishikawa 
et  al.  report  a  TP  cluster  detection  rate  of  85  percent  with  2.5  FPs  per  image  (44,  46). 

Strickland  et  al  (63,  64)  use  40,  50  jtim,  12-bit,  2048x2048  images.  Because  Strick¬ 
land  et  al  used  50  /xm/pixel  images,  this  corresponds  to  a  10x10  cm^  region  on  a  film 
mammogram.  One  of  their  image  pairs  accounts  for  27  percent  of  the  total  clusters  in  the 
database  (104  total  clusters).  Thus,  they  report  results  with  and  without  this  pair.  Over 
all  40  images,  they  achieve  a  TP  detection  rate  of  approximately  84  percent  with  approx¬ 
imately  3  false  clusters  per  image.  (Approximations  are  made  from  the  cluster  detection 
performance  curve  included  in  their  article  (64).)  Without  this  pair,  at  approximately  the 
same  FP  rate,  they  achieve  a  TP  detection  rate  of  approximately  90  percent. 
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Table  1.  Clustered  microcalcification  ROI  detection  results;  NR  indicates  not  reported. 


Investigators 

Database 
(#  images) 

Image  Area 
(cm^) 

Resolution 

(/im/pixel) 

TP  rate 
(%) 

FP  rate 
(#/image) 

Chan  et  al  (overall) 

52 

NR 

100 

92 

1 

Yoshida  et  al. 

39 

8x10 

100 

85 

5 

Nishikawa  et  al 

39 

8x10 

100 

85 

2.5 

Strickland  et  al 

40 

10x10 

50 

84 

3 

2.3  Summary 

As  shown  in  Table  1,  several  candidate  microcalcification  detection  techniques  have 
been  developed  yielding  promising  results;  however,  no  single  method  has  been  developed 
which  is  able  to  reliably  detect  all  microcalcification  clusters.  An  alternative  may  be  to  fuse 
the  results  from  several  complimentary  techniques.  Pattern  recognition  research  at  AFIT 
focuses  on  developing  algorithms  for  use  in  detection  and  classification  of  objects  in  military 
imagery  (17,  26,  37,  51,  52,  59).  Several  of  these  methodologies  have  also  been  applied  to 
the  area  of  breast  cancer  detection  and  diagnosis  (10,  33,  38,  48).  A  new  technique  which 
uses  DoG  filtering  for  detection  of  clustered  micro  calcifications  is  presented  in  the  next, 
chapter. 
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III.  Methodology 


3.1  Introduction 

This  chapter  presents  a  new  technique  for  detecting  and  labeling  clustered  microcal¬ 
cifications  in  digitized  film  mammograms. 

3.2  Algorithm  development 

In  the  literature,  several  CADx  researchers  (9,  44,  46,  68,  64)  develop  and  test  detec¬ 
tion  techniques  on  the  same  database  (not  to  be  confused  with  a  common  database).  This 
research  adopts  the  same  approach.  Although  partitioning  the  database  into  a  training 
and  testing  set  was  considered,  Nishikawa  et  al.  (45)  have  shown  “...the  choice  of  clinical 
cases  used  to  train  and  test  a  computer-aided  diagnosis  (CAD)  scheme  can  affect  the  test 
results  (i.e.  error  rate).”  Nishikawa  et  al.  (45)  also  argue,  “Because  of  the  strong  depen¬ 
dence  of  measured  performance  on  the  testing  database,  it  is  difficult  to  estimate  reliably 
the  accuracy  of  a  CAD  scheme.  Furthermore,  it  is  questionable  to  compare  different  CAD 
schemes  when  different  cases  are  used  for  testing.” 

Although  a  common  CADx  digitized  mammogram  testing  database  has  been  sug¬ 
gested,  this  would  require  standardization  of  pixel  size  and  number  of  gray  levels  used  (45). 
Also,  if  the  common  database  were  enlarged,  all  previously  evaluated  CADx  techniques 
would  need  to  be  reevaluated  if  images  from  new  digitized  imaging  or  film-screen  systems 
were  added  (45).  Furthermore,  Nishikawa  et  al.  argue,  “...because  there  is  no  quantitative 
or  rigorous  means  of  determining  whether  a  particular  common  database  is  representa¬ 
tive  of  the  clinical  population,  it  is  difficult  to  predict  the  true  clinical  impact  of  a  CAD 
scheme.”  Thus,  AFIT  researchers  are  currently  facilitating  a  clinical  study  to  determine 
the  impact  of  this  research  as  well  as  other  AFIT  CADx  schemes  (10,  33,  38,  48). 

3.3  Database 

In  screen-film  mammography,  the  film  is  used  for  detection,  storage,  and  display. 
(15).  Since  full-area  detectors  suitable  for  digital  mammography  are  still  in  the  experi¬ 
mental  stage,  an  alternative  way  of  acquiring  mammograms  in  digital  form  is  to  digitize 
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a  film  mammogram  (67).  This  research  was  performed  using  digitized  film  mammograms 
and  subimages  annotated  by  a  physician  trained  in  radiology  (27). 

Specifically,  the  database  of  56  digitized  (12  bit,  100  //m)  full-breast  (20x10  cm^) 
film  mammograms  contained  63  biopsy-truthed  clustered  microcalcification  ROIs  over  28 
cases.  All  mammograms  included  a  pathology  report  indicating  location  and  diagnosis  of 
biopsied  regions.  The  data  set  included  approximately  2  films  per  case  with  the  following 
distribution:  1  case  with  4  films,  2  cases  with  only  1  film,  and  the  remaining  25  cases  with 
2  films. 

Full-breast  image  optimization  would  have  been  extremely  costly;  approximately  30 
to  60  minutes  per  digitized  film  on  a  Sun  Ultra  or  Sun  Sparc20.  Thus,  52  256x256  and 
9  512x512  subimages  were  extracted  from  these  images  to  facilitate  algorithmic  develop¬ 
ment.  The  relationship  between  full-breast  images,  subimages,  and  ROIs  is  illustrated  by 
example  in  Figure  2.  The  61  subimages  contained  all  63  biopsy-truthed  ROIs,  and  a  physi¬ 
cian  (27)  annotated  abnormalities  identified  in  corresponding  pathology  reports.  These 
subimages  considerably  reduced  algorithmic  development  times,  and  reduced  computer 
memory  requirements.  Additional  information  on  digitization  and  database  management 
is  located  in  Appendix  A. 

S.4  Cluster  detection  criteria 

As  Strickland  et  al.  (64)  note,  no  universally  accepted  criteria  exist  for  the  detection 
of  clustered  microcalcifications.  Thus,  for  each  of  the  subimages  extracted  from  mam¬ 
mograms  with  corresponding  pathology,  a  physician  identifies  the  bounding  box  of  each 
cluster  (as  shown  in  Figure  4a)  and  stores  the  top  left  and  bottom  right  vertices  in  a  truth 
file.  This  technique  for  creating  “ROI  truth”  streamlines  tracking  of  detection  results. 

For  counting  TPs,  a  cluster  is  considered  detected  if  any  cluster  signal  is  found 
within  a  bounding  box  and  the  cluster’s  centroid  is  within  5  mm  from  the  bounding  box. 
A  FP  is  counted  if  all  detections  within  the  cluster  are  not  made  within  a  bounding 
box.  TP  and  FP  detection  rates  are  used  to  determine  receiver  operator  characteristics 
(ROC)  for  performance  measurements  of  radiologic  imaging  systems  (31,  39, 40).  Typically, 
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these  systems  are  optimized  such  that  the  average  TP  rate  is  maximized  while  minimizing 
the  average  FP  rate.  Judy’s  (31)  and  Metz’s  (39,  40)  articles  provide  a  more  complete 
discussion  of  radiologic  ROC  analysis. 

3.5  Noise  reduction  module 

Although  the  main  limitation  in  image  quality  should  be  the  granularity  of  the  film 
emulsion  (67),  noise  is  introduced  from  the  process  of  digitization  (15).  This  noise  may 
be  detected  later  as  a  pseudo  calcification.  There  are  several  noise  reduction  techniques 
available  in  the  literature  (3,  53,  61).  In  this  research,  median  filtering  is  used  because  this 
technique  has  been  shown  to  be  extremely  effective  at  removing  noise  (3).  However,  the 
type  of  median  filter  used  differs  from  the  typical  since  it  is  cross-shaped.  The  cross  shape 
is  formed  by  the  set  of  pixels  which  include  the  center  pixel  and  its  four  nearest  neighbors. 
The  cross  shape  preserves  lines  and  corners  better  than  the  standard  block-shaped  median 
filter  (3,  61).  Additionally,  the  cross  shape  limits  the  possible  substitution  to  the  four 
nearest  neighbors,  reducing  the  potential  for  severe  edge  displacement  (3).  The  subimage 
in  Figure  4a  was  processed  with  this  technique,  and  results  are  shown  in  Figure  4b. 

3.6  Calcification  enhancement  module:  DoG  filter 

After  noise  has  been  reduced,  the  image  is  DoG  filtered.  This  section  discusses  DoG 
filtering  by  example.  First,  the  physiological  basis  for  and  construction  of  a  DoG  filter  is 
identified.  Next,  DoG  spatial  kernels  are  generated  based  on  potential  target  sizes.  These 
filters  are  then  individually  applied  to  the  noise  reduced  subimage  shown  in  Figure  4b. 

3.6.1  Physiological  basis.  In  chapter  two  of  Vision  (36),  Marr  reviews  require¬ 
ments  for  and  efficient  detection  of  intensity  changes,  identifies  the  physiological  basis  for 
the  DoG,  and  (indirectly)  provides  an  engineering  rule  of  thumb  for  construction  of  DoG 
filters.  To  detect  intensity  changes  efficiently,  Marr  argues  one  should  search  for  a  filter 
which  is  a  differential  operator  and  is  tunable  to  act  at  any  desired  scale  (36).  Raff  (49) 
provides  further  insight  into  Marr’s  argument.  In  short,  two  basic  ideas  underlie  detection 
of  intensity  changes  in  images.  First,  intensity  changes  occur  at  different  scales  in  images. 
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Figure  4.  256  x  256  subimage:  (a)  ROI  bounding  box  contains  biopsy-truthed  malignant 

cluster  of  microcalcifications;  (b)  noise-reduced  subimage 
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and  second,  sudden  changes  in  intensities  result  in  peaks  of  first  derivatives  and  hence 
“zero  crossings”  in  the  second  derivatives. 

Marr  argues  the  Laplacian  of  Gaussian  (LoG)  is  the  optimum  filter  which  satisfies 
these  conditions  (36).  There  are  two  basic  steps  to  LoG  filtering.  First,  use  a  Gaussian 
to  effectively  wipe  out  all  structure  at  scales  smaller  than  the  space  constant  (standard 
deviation)  of  the  Gaussian.  Second,  use  a  Laplacian  as  a  differential  operator  to  detect 
intensity  changes  in  the  blurred  image  (36).  As  Marr  notes,  the  benefit  of  using  a  Gaussian 
is  that  it  is  smooth  and  localized  in  both  the  space  and  frequency  domains.  Additionally, 
the  benefits  of  using  the  Laplacian  are  that  it  is  an  orientation-independent  operator,  and 
it  reduces  the  extra  computational  burdens  associated  with  using  second-order  directional 
derivative  operators  (36).  Next,  Marr  identifies  the  biological  implications  of  the  LoG 
filter. 

Neurophysiological  experiments  provide  evidence  the  human  visual  pathway  includes 
a  set  of  “channels”  that  are  orientation  and  spatial  frequency  selective  (36,  49,  61).  Essen¬ 
tially,  at  each  point  in  the  visual  field,  there  are  four  size- tuned  filters  or  masks  analyzing 
the  image  (36).  The  operation  of  these  spatial  receptive  fields  (in  the  form  of  retinal  gan¬ 
glion  cells)  can  be  approximated  closely  by  a  DoG  (36,  49,  61).  Marr  argues  LoG  filters 
form  the  basis  for  these  psychophysically  determined  channels  and  the  LoG  operator  can 
be  approximated  closely  by  a  DoG  when  the  ratio  of  the  DoG  space  constants  is  1:1.6  (36). 

3.6.2  DoG  ■parameters.  The  2-D  Gaussian  smoothing  operator  is  given  as 

—  (  )  /  \ 

=  ce  (1) 

where  c  normalizes  the  sum  of  mask  elements  to  unity,  x  and  y  are  horizontal  and  vertical 
indices,  and  a  is  the  space  constant  (61).  Using  Equation  1,  the  difference  of  two  Gaussians 
with  different  a  yields 

DoG{x,y)  =  Cie  —  C2e  (2) 
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For  a  target  of  size  (average  width)  t  pixels,  use 


(3) 


and,  from  Marr’s  ratio  (36),  use 


(4) 


Figure  5.  DoG  kernels  with  their  corresponding  frequency  response;  (a)  t  =  1.0,  (b) 
t  =  2.0,  and  (c)  t  =  3.0.  Units  of  t  are  pixels,  where  1  pixel  =  100  firsi. 


3.6.3  DoG-filtered  examples.  Since  microcalcifications  typically  range  from  100 
to  300  /im  in  diameter  (8,  56),  potential  target  sizes  for  the  100  /im  digitized  mammo¬ 
grams  correspond  to  1  to  3  pixels.  Using  t  =  {1.0, 2.0, 3.0},  we  find  a2  =  {0.5, 1.0, 1.5} 
and  CTi  =  {0.3125,0.6250,0.9375},  respectively.  DoG  filters  were  constructed  using  these 
paired  parameters  as  shown  in  Figure  5.  These  kernels  were  individually  used  as  spatial 
filters  on  the  noise-reduced  subimage  shown  in  Figure  4b  with  corresponding  results  shown 
in  Figure  6.  Alternatively,  an  “optimized”  DoG  kernel  is  shown  in  Figure  7a  with  corre- 
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(a)  (b)  (c) 


Figure  6.  Results  of  applying  three  DoG  kernels  to  a  noise-reduced  subimage;  (a)  t  =  1.0, 
(b)  t  =  2.0,  and  (c)  t  =  3.0.  Units  of  t  are  pixels,  where  1  pixel  =  100  /xm.  As 
taxgetsize  increases,  larger  image  structures  are  emphasized  as  expected. 


Urgetalze-1 .20 
DoG  filter  coefficients 


Figure  7.  Results  of  DoG  filtering  noise  reduced  subimage  using  “optimized”  targetsize; 
(a)  DoG  kernel,  t  =  1.1995  pixels,  and  (b)  DoG  filtered  subimage 
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spending  filtered  subimage  results  shown  in  Figure  7b.  (Optimization  of  t  is  discussed  in 
Section  3.10.) 

S.7  Calcification  detection  module 

Once  an  image  has  been  preprocessed  to  enhance  differences  between  targets  and 
background,  the  image  is  more  easily  segmented.  The  DoG-filtered  subimages  (Figure  6  and 
Figure  7b)  contain  obvious  differences  in  gray  levels  between  potential  micro  calcifications 
and  background;  however,  the  optimum  method  for  segmentation  is  not  so  clear.  Although 
micro  calcifications  tend  to  be  among  the  brightest  objects  in  DoG-filtered  subimages,  they 
may  exist  within  regions  of  high  average  gray  levels  and  thus  prove  difficult  to  reliably 
segment.  A  method  that  generally  addresses  these  concerns  uses  pair-wise  pixel  “anding” 
of  the  results  of  global  histogram  and  locally  adaptive  thresholding. 

If  targets  tend  to  exist  within  an  image’s  higher  gray  levels,  then  the  global  threshold 
may  be  approximated  by  finding  the  level  which  segments  a  preselected  percentage  of  the 
corresponding  higher  levels  in  the  image  histogram.  Locally  adaptive  thresholding  may  be 
implemented  by  varying  the  high  and  low  thresholds  based  on  the  local  pixel  value  mean 
and  standard  deviation  (43). 

According  to  Gonzales,  “the  histogram  of  a  digital  image  with  gray  levels  in  the 
range  [0,  L-1]  is  given  as 

p{rk)  =  —  (5) 

n 

where  r*.  is  the  fcth  gray  level,  rijs,  is  the  number  of  pixels  in  the  image  with  that  gray  level, 
n  is  the  total  number  of  pixels  in  the  image,  and  k  =  0, 1,2,...,L  —  1”  (24).  After  the 
histogram  has  been  computed,  the  gray  level,  g,  used  to  segment  (threshold)  a  preselected 
upper  fraction  of  the  histogram,  /,  is  found  using 

/  =  1  -  Ylpi'^k)  (6) 

k-o 

where  0  <  g  <  L  —  1, 
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The  locally  adaptive  thresholds,  tio  and  tku  found  using 


tio  =  kio(r{x,  y)  +  m{x,y)  (7) 

and 

thi  =  khi'y{x,y)  +  rn{x,y)  (8) 

where  kio  and  ku  are  used  to  preselect  the  multiple  of  a{x,y),  the  local  standard  deviation, 
and  m{x,y)  is  the  local  mean  of  the  iV  x  iV  neighborhood  centered  on  the  pixel  at  (x,j/) 
(43).  (Optimization  of  /,  ku^  khi,  and  N  is  discussed  in  Section  3.10.) 

Once  the  globally  and  locally  thresholded  images  are  generated,  they  are  combined 
into  an  initial  calcification  detection  mask,  Mj,  using 

Mi{x,y)  =  Tg{x,y)r\Tio{x,y)r\Thi{x,y)  (9) 

where  Tg{x,y)  is  the  globally  thresholded  image,  and  Tio{x,y)  and  Thi{x,y)  are  the  images 
locally  thresholded  using  kia  and  kh{  respectively.  For  example,  this  method  was  performed 
on  the  DoG  filtered  subimage  in  Figure  7b  with  results  shown  in  Figure  8a.  To  facilitate 
clustering,  the  centroids  of  individual  detections  are  used  to  generate  the  final  calcification 
detection  mask,  M/,  shown  in  Figure  8b. 

When  working  with  full-breast  images,  however,  detections  may  occur  on  film  labels 
which  are  obvoiusly  not  ROIs.  Thus,  to  ensure  detections  are  within  the  breast,  a  breast 
segmentation  utility  was  developed  to  segment  the  breast  from  the  image.  This  utility 
consists  of  five  steps.  First,  a  variance  image  is  computed  from  a  subsampled  version  of 
the  full  breast  image.  Second,  the  original  image  and  the  variance  image  are  rescaled 
and  added.  Next,  Floyd-Steinberg  image  dithering  is  used  to  reduce  the  grayscale  image 
to  two  gray  levels.  Only  the  largest  contiguous  region  is  set  to  one  assuming  it  always 
corresponds  to  the  breast.  All  other  regions  are  set  to  zero.  Finally,  the  resulting  binary 
image  is  upsampled  to  generate  a  full-image  breast  mask.  This  utility  returns  a  binary 
breast  mask  with  ones  corresponding  to  the  region  of  film  containing  the  breast. 
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Thus,  to  ensure  detections  are  within  the  breast,  the  centroids  of  individual  locations 
are  compared  with  the  breast  mask.  If  the  centroids  are  found  within  the  breast  mask,  the 
centroids  are  then  passed  to  the  cluster  detection  module. 

3.8  Cluster-detection  module 

Recall,  “...malignant  calcifications  tend  to  occur  in  tight  clusters  of  1  cm  diameter 
or  less  ”  (56);  however,  calcifications  can  be  extensive  (56).  Thus,  the  cluster  detection 
module  identifies  clusters  based  on  the  clustering  hypothesis  proposed  by  Chan,  et  al.  (9). 
Specifically,  a  suspicious  ROI  (cluster)  is  fJ-Csmin  or  more  detected  signals  linked  by  an 
Euclidean  nearest  neighbor  distance,  dnn-  (Optimization  of  i^Csmin  ^nd  is  discussed 
in  Section  3.10.) 

3.9  Labeling  computer- detected  ROIs 

Once  detections  have  been  indexed  according  to  their  corresponding  ROIs,  an  ellipse 
may  be  used  to  identify  computer-detected  ROIs,  as  shown  in  Figure  9a.  The  shape  of  the 
ellipse  may  be  generated  based  on  the  distribution  of  cluster  signals.  Figure  9b  provides  an 
example  of  the  final  output  with  corresponding  truth  for  the  subimage  shown  in  Figure  4a. 
Figure  10  shows  the  corresponding  full  image  result  which  may  be  provided  in  hardcopy 
form  to  the  radiologist.  Figure  11  illustrates  overall  detection  methodology. 
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Figure  9.  Cluster  detection  results  for  t  =  1.1995,/  =  0.0009,  =  4.1226, /c/ii  = 

7.000,  iV  =  69^tiCsmin  =  3,  and  d^n  =  51.0031;  (a)  computer-detected  ROI 
with  corresponding  microcalcification  detections,  and  (b)  computer-detected 
ROI  for  radiologist  consumption.  Since  the  “ground  truth”  for  the  cluster  was 
available,  the  corresponding  bounding  box  is  shown  and  detection  results  are 
reported. 
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Figure  10.  Full  image  (20x10  cm^)  cluster  detection  results  for  t  =  1.1995,/  = 
0.0009,  ku  =  4.1226,  hi  =  7.000,  N  =  69,  nCsmin  =  3,  and  d„„  =  51.0031; 
Since  the  “ground  truth”  for  the  cluster  was  available,  the  corresponding 
bounding  box  is  shown.  The  upper  ROI  not  corresponding  to  the  truth  box 
contained  a  true  cluster  of  microcalcifications,  however  the  region  was  not 
biopsied.  The  lowest  ROI  was  due  to  scratches  on  the  film. 
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Figure  11.  Cluster  Detection  Methodology:  t  controls  the  DoG  filter  used  to  enhance 
calcifications;  /,  fc;„,  khi,  and  N  control  the  detection  stage  of  the  system; 
IJ’Csmin  and  control  the  clustering  module. 
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3 AO  Cluster  detection  optimization 


As  shown  in  Figure  11,  the  following  seven  parameters  are  used  to  control  clustered 
microcalcification  detection: 

•  t,  microcalcification  target  size  (average  width)  in  pixels  (1  pixel  =  100  /xm);  i.e., 
t  =  {1.0, 1.15, 2.13, . . .},  where  t  >  0 

•  /,  upper  fraction  in  intensity  oftheDoG-filtered  image  histogram;  i.e.,/  =  {0.001,0.02,0.05,...}, 
where  0  <  /  <  1 

•  kio  and  khi,  multiples  of  cr{x,  y),  the  standard  deviation  in  brightness  computed  from 

neighborhood  centered  on  DoG  filtered  pixel  at  (x,y);  i.e.,  =  {3.6,3.75,4.17, . . .}, 

and  khi  =  {6.3,7.0,8.25,...}  where  0  < 

•  N,  block  size  of  neighborhood  centered  on  DoG-filtered  pixel  at  (x,y);  i.e.,  N  = 

{51, 53, 67, . . .},  where  JV  >  0,  and  iV  €  Odd 

•  fiCsmim  minimum  number  of  detected  microcalcifications  linked  by  to  trigger  a 
ROI  detection;  i.e.,  fiCsmin  =  {2,3,4,5,. . .},  where  /J-Csmin  >  2 

•  d„„,  Euclidean  nearest  neighbor  distance  (in  pixels)  between  detected  microcalcifica¬ 
tions;  i.e.,  dnn  -  {50.23,51.5,60.4, . . .},  where  >  0 

Recall,  t  controls  the  DoG  filter  used  to  enhance  calcifications.  Additionally,  /, 
kio,  khi,  and  N  control  the  detection  stage  of  the  system.  Finally,  dnn  and  piCsmin  control 
the  clustering  module.  Although  empirical  results  may  be  obtained  by  manually  modifying 
these  parameters,  an  alternative,  more  efficient  method  is  to  use  an  optimization  algorithm. 


3.10.1  Function  optimization.  A  function  optimization  problem  may  be  described 
as  follows:  given  some  finite  domain,  D,  a  particular  setting  of  cluster  detection  parameters, 
X  =  {t,f,kio,khi,N,iJ.Csmin,dnn}  where  x  G  D,  and  an  objective  function  fotj  :  D  3?, 
find  the  x  in  D  that  maximizes  or  minimizes  fobj  (61).  Optimization  problems  require  an 
objective  function  which  reflects  the  “goodness”  of  a  solution  (61). 


fobj{-x)  = 


-FP{x)  ,  rP(x)  >  TPmi. 
FPpenaity  ,  Otherwise 


(10) 
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Recall,  radiologic  imaging  systems  may  be  optimized  to  maximize  the  TP  rate  while 
minimizing  the  FP  rate.  This  objective  may  be  recast  into  the  functional  form  shown  in 
Equation  10,  where  maximizationis  the  goal.  (The  genetic  algorithm  technique  used  in  this 
research  seeks  maximization  of  an  objective/fitness  function.)  For  a  particular  setting  of 
cluster  detection  parameters,  if  the  minimum  acceptable  TP  rate,  TPmin,  is  exceeded,  the 
objective  function  returns  the  negative  of  the  FP  rate.  Otherwise,  if  the  TP  rate  falls  below 
TPmin,  the  objective  function  returns  a  constant  value,  FPpenaity',  i-e-,  FPpenaity  =  -10 
because  — F’P(x)  values  observed  in  this  research  typically  range  from  approximately  —.27 
to  -1.0.  Although  Richardson  et  al.  (50)  argue  “...infeasible  solutions  should  provide 
information  and  not  just  be  thrown  away,”  they  acknowledge  “...many  researchers  believe 
penalty  functions  should  be  harsh”  so  optimization  methods  “...will  avoid  the  forbidden 
spaces.”  Furthermore,  since  Richardson  et  al  (50)  support  their  argument  with  only  two 
contrived  optimization  problems,  a  harsh  FPpenaity  is  used. 

Conventional  approaches  to  optimization  use  calculus-based  methods  (61).  However, 
these  approaches  are  generally  limited  to  objective  functions  that  possess  “nice”  properties 
such  as  continuity,  differentiability,  etc.  (29)  Their  main  limitation  is  their  local  behavior; 
the  search  can  easily  end  in  a  local  maximum,  and  the  global  maximum  can  be  missed 
(61).  Thus,  several  techniques  are  used  to  improve  the  probability  of  finding  the  global 
maximum:  start  at  several  points  in  the  search  space,  use  dynamic  programming,  apply 
random  searches,  etc.  (61).  One  relatively  new  and  widely  applicable  approach  is  to  use 
genetic  algorithms  augmented  by  one  of  these  standard  methods  (28,  61). 

3.10.2  Genetic  algorithms.  Genetic  algorithms  (GAs)  have  been  successfully 
applied  to  many  diverse  and  difficult  optimization  problems  (29)  such  as  jet  engine  tur¬ 
bine  design  (28),  tracking  of  criminal  suspects  through  face  composites  (6),  and  binary 
phase-only  filter  design  for  pattern  recognition  (7).  GAs  “...are  at  their  best  when  explor¬ 
ing  complex  landscapes  [solution  spaces]  to  find  regions  of  enhanced  opportunity”  (28). 
However,  as  Holland  (inventor  of  the  GA)  (28)  notes,  “...if  a  partial  solution  can  be  im¬ 
proved  further  by  making  small  changes  in  a  few  variables,  it  is  best  to  augment  the  genetic 
algorithm  with  other,  more  standard  methods.”  Thus,  this  research  uses  a  Matlab  imple- 
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mentation  of  a  GA  developed  by  Houck  et  al  (29)  to  find  promising  parameter  settings, 
then  augments  the  GA  with  the  simplex  method  (62)  included  with  Matlab. 

GAs  search  the  solution  space  to  maximize  a  fitness  (objective)  function  using  simu¬ 
lated  evolutionary  operators;  i.e.,  mutation,  sexual  recombination,  etc.  (28,  29,  61).  GAs 
manipulate  and  maintain  a  population  of  solutions  while  searching  for  6e«er  solutions  (29). 
The  fittest  individuals  of  any  population  tend  to  reproduce  more  and  survive  longer.  This 
increases  the  average  fitness  in  successive  generations  (29,  61).  However,  as  Houck  et  al. 
(29)  note,  inferior  individuals  can  survive  and  also  reproduce.  A  more  complete  discussion 
of  GAs  and  related  topics  can  be  found  in  Michalewicz’s  book  (41).  Houck  (29)  summarizes 
the  GA  as  follows: 

1.  Supply  [an  initial]  population  Pq  of  [n]  individuals  and  respective  [fitness]  function 
values 

2.  z  ^  1 

3.  P'i  •t—  selection.function{Pi  —  1) 

4.  Pi  *-  reproduction-function{Pl) 

5.  z  z  -k  1 

6.  Repeat  step  3  until  termination  [termination  criteria  defined  later] 

7.  Print  out  best  solution  found 

GA  use  requires  determination  of  several  issues:  objective  function  design,  param¬ 
eter  set  representation,  population  initialization,  choice  of  selection  function,  choice  of 
genetic  operators  (reproduction  mechanisms)  for  simulated  evolution,  and  identification 
of  termination  criteria  (29,  61).  Although  one  may  argue  the  GA  itself  appears  to  re¬ 
quire  optimization  of  GA  parameter  settings,  several  researchers  have  closely  examined 
this  “problem”  (4,  11,  22,  23,  30,  50,  54).  As  there  is  no  standard  agreement  on  optimum 
GA  parameter  settings,  this  research  adopts  the  following  general  “rules  of  thumb”  culled 
from  several  investigations  in  the  literature: 
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•  As  Sonka  et  al.  (61)  note,  “No  optimization  problem  can  guarantee  finding  a  good 
solution  to  the  problem  if  the  objective  function  does  not  reflect  the  goodness  of 
the  solution.  Therefore,  the  design  of  the  objective  function  is  a  key  factor  in  the 
performance  of  any  optimization  algorithm.”  Recall  Equation  10:  if  an  acceptable 
TP  rate  is  reached,  then  minimization  of  the  FP  rate  is  desired. 

•  For  parameter  representation,  Janikow  and  Michalewicz  (30,  41)  showed  the  real¬ 
valued  GA  is  an  order  of  magnitude  more  efficient  in  CPU  time  than  the  binary  GA 
and  provides  higher  precision  with  more  consistent  results  across  replications  (29). 
Thus,  a  floating-point  representation  is  used  in  this  research. 

•  Because  GAs  can  iteratively  improve  existing  solutions  (29),  several  GA  users  seed 
the  initial  population  with  some  members  known  beforehand  to  be  in  an  interesting 
part  of  the  search  space  (4,  22).  This  technique  is  used  because  preliminary  cluster 
detection  parameter  settings  appear  to  be  in  an  interesting  region  of  the  search 
space.  Also,  while  the  minimum  number  of  members  in  a  GA  population  is  three 
(22),  Schaffer  et  al.  (54)  argue  good  online  performance  can  be  expected  with  20  to 
30  members.  Because  Schaffer  et  al.  support  their  suggestion  using  a  large  problem 
suite  over  8400  genetic  searches  (lOK  evaluations/search)  (54),  and  each  objective 
function  evaluation  is  computationally  expensive  (approximately  30  to  60  minutes 
on  a  Sun  Ultra  or  Sun  Sparc20)  20  members  are  used  in  this  research. 

•  Several  schemes  exist  for  the  (typically)  probabilistic  selection  process  used  to  identify 
candidates  for  reproduction:  roulette  wheel  or  proportional  selection,  scaling  tech¬ 
niques,  tournament,  elitist  models,  and  ranking  methods  (29,  41,  61).  According  to 
Grefenstette  and  Baker  (25),  although  several  selection  algorithms  have  been  studied, 
many  “...have  not  enjoyed  the  same  level  of  theoretical  analysis  as  the  proportional 
selection  algorithm.”  The  authors  argue  a  promising  alternative  to  proportional 
selection  is  ranking  because  “the  resulting  algorithm  is  less  prone  to  premature  con¬ 
vergence  cause  by  individuals  that  are  far  above  average.”  The  basic  idea  of  ranking 
is  to  select  strings  (solutions)  for  the  mating  pool  based  on  the  relative  fitness  be¬ 
tween  strings.  This  research  uses  the  ranking  method  included  in  the  Matlab  GA, 
normalized  geometric  ranking  (29). 
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•  Two  basic  types  of  genetic  operators  are  commonly  used  in  GAs:  crossover  and 
mutation  (28,  54,  61).  According  to  Sonka  et  al  (61),  “...the  basic  idea  [of  crossover] 
is  to  randomly  mate  the  . . .  strings  [solutions]  ”  identified  by  the  selection  process  by 
“...randomly  choosing  a  position  for  the  border  of  each  pair  of  strings,”  and  producing 
“...new  strings  by  swapping  all  characters  between  the  beginning  of  the  string  pairs 
and  the  border  position.”  Sonka  et  al.  (61)  note  the  mutation  operator’s  role  is 
to  “...randomly  change  one  character  of  some  string  of  the  population  from  time  to 
time.”  Mutation  is  used  because  “...some  local  configuration  of  characters  in  strings 
of  the  population  can  be  totally  lost”  due  to  the  selection  and  crossover  operations 
(61).  As  there  is  no  general  agreement  on  specific  crossover  and  mutation  schemes, 
this  research  adopts  the  default  crossover  and  mutation  schemes  included  with  the 
Matlab  GA:  arithmetic  crossover  and  nonuniform  mutation  (29). 

•  Finally,  the  termination  criteria  used  in  this  research  adopts  the  commonly  used 
technique  of  running  each  iteration  until  a  specified  number  of  generations  are  done 
(29,  54);  here,  ten  generations  are  used  due  to  the  high  computational  cost  men¬ 
tioned  earlier.  Although  termination  due  to  loss  of  population  diversity  and/ or  lack 
of  improvement  is  efficient  when  crossover  is  the  primary  source  of  variation  in  a  pop¬ 
ulation  (4,  22,  29),  homogeneous  populations  can  be  succeeded  by  populations  with 
better  (higher)  fitness  when  using  mutation  (4).  Additionally,  since  Goldberg  (22) 
argues  restarting  populations  that  may  have  converged  proves  useful,  this  research 
rims  several  iterations  of  the  GA  until  a  consistent  lack  of  improvement  is  recognized. 

3.10.3  Simplex  method.  As  mentioned  earlier,  once  potentially  optimum  solu¬ 
tions  are  found  by  the  GA,  the  most  fit  GA  solution  may  be  further  optimized  using  the 
simplex  method  developed  by  Dantzig  (28,  62).  The  Matlab  implementation  of  the  sim¬ 
plex  algorithm  requires  a  basic  (initial)  feasible  solution,  Xq  (found  by  the  GA),  and  a 
cost  function  (Equation  11)  (5,  62).  In  this  research,  the  simplex  method  is  used  to  locate 
other  feasible  solutions  (members  of  the  “feasible  set”)  while  minimizing  the  cost  func¬ 
tion  until  the  optimal  solution  is  obtained  (5).  It  should  be  noted,  since  the  GA  provides 
potentially  optimum  solutions  not  necessarily  near  or  equal  to  the  global  optimum,  the 
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simplex  method  is  not  guaranteed  (nor  expected)  to  locate  the  global  optimum.  Rather, 
the  simplex  method  is  used  to  further  optimize  the  most  fit  solution  found  by  the  GA.  A 
cursory  overview  of  the  simplex  method  follows. 


First,  several  terms  must  be  defined.  According  to  Bronson  (5),  “A  mathematical 
program  is  an  optimization  problem  in  which  the  objective  and  constraints  are  given  as 
mathematical  functions  and  functional  relationships.”  If  these  relationships  are  linear  in 
each  of  their  arguments  a  mathematical  program  is  linear  (5).  The  simplex  method  was 
developed  to  solve  linear  programs  (62).  Although  the  cost  function  and  its  constraints 
(Equation  11)  may  not  be  strictly  linear,  the  simplex  method  is  effective  at  further  opti¬ 
mizing  the  most  fit  solutions  found  by  the  GA. 

The  cost  function  (Equation  11)  brings  to  the  problem  a  family  of  parallel  planes 
(62).  The  solution  with  zero  cost  is  the  plane  that  goes  through  the  origin;  other  planes 
give  all  other  possible  costs  (62).  These  planes  sweep  out  the  entire  n-dimensional  solution 
space,  and  the  optimal  solution  occurs  where  the  planes  first  touch  the  feasible  set  (62). 

Thus,  as  Strang  (62)  notes,  there  are  basically  the  two  phases  in  the  simplex  method: 


1.  Locate  a  corner  of  the  feasible  set,  where  a  corner  is  the  meeting  point  of  n  different 
planes. 

2.  Within  the  feasible  set,  go  from  corner  to  corner  along  an  edge  guaranteed  to  decrease 
the  cost.  Terminate  when  a  corner  is  reached  from  which  all  edges  increase  the  cost. 

Finally,  while  the  total  number  of  simplex  steps  needed  to  find  the  optimum  solu¬ 
tion  is  impossible  to  answer,  Strang  (62)  reports  the  simplex  method  has  been  shown  to 
operate  in  an  average  time  which  is  polynomial  (as  opposed  to  the  undesirable  class  of  non¬ 
polynomial,  or  NP).  Because  each  cost/objective  function  evaluation  is  computationally 
expensive  (approximately  30  to  45  minutes  on  Sun  Ultra  or  Sun  Sparc20),  this  research 
uses  300  simplex  steps.  A  more  complete  discussion  of  the  simplex  method  may  be  found 
in  Strang  (62). 
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3.11  Summary 

This  chapter  presented  a  new  technique  to  detect  clustered  microcalcifications  in 
digitized  film  mammograms.  First,  the  image  is  preprocessed  to  reduce  digitizaion  noise. 
Next,  a  DoG  filter  is  used  to  enhance  microcalcifications.  The  filtered  image  is  thresholded, 
and  detections  are  passed  to  the  cluster  detection  module.  Finally,  computer-detected 
ROIs  are  overlaid  on  the  original  image  for  radiological  analysis.  To  improve  detection 
rates,  globally  optimize  detection  parameters  using  a  genetic  algorithm  (GA),  then  locally 
optimize  using  the  simplex  method. 


36 


IV.  Results 


This  chapter  presents  the  results  of  the  clustered  micro  calcification  detection  tech¬ 
nique. 

4.1  Optimization 

As  mentioned  earlier,  full-breast  image  optimization  of  detection  parameters  would 
have  been  computationally  expensive.  Thus,  subimages  were  used  during  the  optimization 
stage.  All  subimages  had  at  least  one  biopsy-truthed  ROI.  If  at  least  one  subimage  ROI 
was  detected,  a  TP  of  1  was  recorded.  All  FPs  within  each  subimage  were  recorded  as  well. 
The  goal  of  the  optimization  stage  was  to  reduce  the  subimage  FP  rate  while  maintaining 
the  minimum  acceptable  TP  rate,  TPmin  =  84.5  percent. 

4.1.1  TPmin-  Since  the  optimization  method  described  in  Chapter  III  requires 
a  minimum  acceptable  TP  rate  to  minimize  the  FP  rate,  an  optimized  ROC  analysis  was 
performed  using  TPmin  =  84.5,82.5,80.0,  and  77.0  percent.  Although  higher  subimage 
TPmin  goals  were  easily  attainable,  they  were  accompanied  by  higher  subimage  FP  rates 
as  expected;  approximately  0.8  to  2.6  FPs  per  subimage. 

Recall,  image  area  is  an  important  factor  in  increasing  the  FP  rate  since  the  larger  the 
test  area,  “...the  greater  the  chance  for  producing  false  positives  in  general”  (18).  Thus, 
because  most  of  the  subimages  were  approximately  1/32  the  area  of  a  database  image, 
subimage  detection  rates  of  0.8  to  1.3  may  lead  to  approximately  25.6  to  41.6  FPs  per 
full-breast  image.  Continuing  with  the  analysis,  if  a  maximum  of  approximately  5  FPs  per 
full-breast  image  is  desired,  approximately  0.1562  FPs  per  subimage  may  be  considered 
the  goal. 

4.2  Subimage  ROC  results 

To  provide  insight  from  a  ROC  perspective  into  how  subimage  TP  and  FP  detection 
rates  vary  with  the  choice  of  TPmini  the  GA  and  simplex  methods  were  used  to  search 
the  feasible  set  of  parameter  settings.  Additionally,  the  case  TP  rate  was  tracked  to 
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Table  2.  TPmin  vs.  Optimized  subimage  parameter  settings 


TPmin 

t 

/ 

kio 

^hi 

N 

dnn 

l^Cs>fYiin 

84.5 

1.1995 

0.0009 

4.1226 

7.0000 

69 

51.3019 

3 

82.5 

1.1995 

0.0009 

4.1226 

7.0000 

69 

51.3420 

3 

80.0 

1.2286 

0.0009 

7.0000 

69 

51.3019 

3 

77.0 

1.2276 

0.0007 

3 

Table  3.  Optimized  subimage  ROC  results 


TPmin 

Case  TP 

Subimage  TP 

Subimage  FP 

96.43 

85.25 

0.2295 

85.25 

81.97 

77.0 

92.86 

81.97 

0.1803 

demonstrate  optimization  effects.  Optimized  subimage  parameter  settings  found  during 
this  search  are  summarized  in  Table  2,  with  corresponding  TP  and  FP  rates  in  Table  3. 
Recall,  t  controls  the  DoG  filter  used  to  enhance  calcifications.  Additionally,  /,  kio,  khi, 
and  N  control  the  detection  stage  of  the  system.  Finally,  and  i^Csmin  control  the 
clustering  module. 

As  TPmin  was  lowered,  the  subimage  FP  rates  fell  as  expected,  however  the  corre¬ 
sponding  drop  in  subimage  TP  rates  was  accompanied  by  a  drop  in  the  case  TP  rate  from 
96.43  to  92.86  percent.  Also,  as  mentioned  earlier,  although  higher  subimage  TPmin  goals 
were  attainable,  they  were  accompanied  by  higher  subimage  FP  rates  as  expected.  Thus, 
optimized  parameter  settings  found  using  TPmin  =  84.5  percent  were  chosen  for  full  digi¬ 
tized  mammogram  image  evaluation.  These  settings  corresponded  to  a  TP  case  detection 
rate  of  96.43  percent  and  a  subimage  detection  rate  of  85.25  percent  with  0.2295  FPs  per 
subimage. 

4.3  Database  results 

4.3.1  Physician  Evaluation  of  Detection  Scheme.  Ultimately,  full  image  results 
were  evaluated  as  shown  in  Table  4  using  the  following  optimized  parameter  settings: 
t  =  1.1995,/  =  0.0009,  fcjo  =  4.1226,  few  =  7.000,  AT  =  69,  and  laCsmin  =  3.  To  ensure  full 
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image  detection  results  were  accurate,  a  physician  trained  in  radiology  (27)  individually 
verified  each  ROI  detection  on  the  original  film  as  a  TP  or  FP.  (Recall,  the  physician  only 
truthed  subimages,  not  full  images.)  Again,  if  at  least  one  image  ROI  was  detected,  a 
TP  of  1  was  recorded.  Also,  since  the  detection  algorithm  was  evaluated  on  full-breast 
images,  several  additional  non-biopsied  microcalcification  ROIs  were  detected.  Thus,  an 
additional  category  of  “Xtra  ROIs”  was  used  to  track  additional  physician  verified  ROIs. 
Recall,  Figure  10  shows  an  example  of  an  “Xtra  ROI.”  All  FPs  within  each  subimage  were 
recorded  as  well. 

To  provide  insight  into  the  sources  of  FPs,  the  physician  (27)  used  the  following 
two  additional  categories  to  describe  noteworthy  FPs:  Scratch  and  BB  induced  ROIs.  As 
mentioned  earlier,  scratches  and  pick-off  of  film  emulsion  are  a  known  source  of  pseudocal¬ 
cifications  in  film  mammograms  (20).  The  term  “scratch  induced  ROF  is  used  for  ROIs 
due  to  scratches  and/or  film  emulsion  pick-off.  Recall,  Figure  10  shows  an  example  of 
a  scratch  induced  ROI.  BBs  are  typically  used  as  a  marker  on  mammograms  to  identify 
palpable  lesions  in  the  breast.  It  should  be  noted  that  no  TPs  were  due  to  BBs.  Also, 
as  shown  in  Table  4,  only  9  of  56  images  had  double-digit  FP  rates.  The  high  FP  rate 
on  these  images  may  be  attributed  to  parenchymal  edge  effects;  the  functioning  structure 
of  these  breasts  (ducts,  lobules,  etc.)  introduced  a  large  number  of  thin,  sharp  densities 
which  appeared  calcification-like  to  the  detection  scheme. 

4.3.2  Analysis:  Detection  Scheme  Results.  The  detection  technique  demon¬ 
strated  the  following  results  (Table  4)  over  56  full-breast  digitized  film  mammograms: 

•  96.4%  (27/28)  TP  case  detection  rate 

•  85.7%  (54/63)  TP  ROI  detection  rate 

•  85.7%  (48/56)  TP  image  detection  rate 

•  5.75  FPs  per  full-breast  image  (20x10  cm,  12  bit,  100  /xm) 

•  0.68  scratch  induced  FPs  per  image 

•  3  BB  induced  false  ROIs 

•  2.1964  “Xtra  ROIs”  per  image 


39 


Table  4.  Results  using  t  =  1.1995,/  =  0.0009,  fcjo  =  4.1226,  =  7.000,  iV  =  69,  and 

=  3:  a  TPca.e  rate  of  96.43%,  and  TProi  and  TPi„,age  rates  of  85.7% 
with  5.75  FPs  per  full-breast  image.  Note:  Only  9  images  had  double  digit  FPs 
mostly  due  to  parenchymal  edge  effects. 
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Table  5.  Clustered  microcalcification  ROI  detection  results;  NR  indicates  not  reported. 


Investigators 

Number 

of 

images 

Image 

Area 

(cm^*) 

Resolution 

(;tim/pixel) 

TP 

(%) 

FPs 

(#/image) 

Area  ad¬ 
justed 

FPs 

(#/image) 

Chan  et  al  (overall) 

52 

NR 

100 

92 

1 

NR 

Yoshida  et  al 

39 

8x10 

100 

85 

5 

12.5 

Nishikawa  et  al 

39 

8x10 

100 

85 

2.5 

6.25 

Strickland  et  al 

40 

10x10 

50 

84 

3 

6 

Ochoa  et  al 

56 

20x10 

100 

85.7 

5.75 

5.75 

As  mentioned  earlier,  a  direct  comparison  of  detection  results  against  other  CADx 
detection  schemes  (9,  44,  46,  64,  68)  cannot  be  done  since  CADx  researchers  use  different 
databases.  However,  as  shown  in  Table  5,  the  85.7  percent  ROI  TP  detection  rate  found 
using  this  detection  scheme  is  comparable  with  reported  ROI  TP  detection  rates  of  84  to 
93  percent  (9,  44,  46,  64,  68).  Although  the  performance  of  a  CAD  scheme  is  difficult 
to  reliably  estimate  (45),  the  number  of  images  used  to  determine  TP  and  FP  rates  was 
larger  than  those  used  by  Chan  et  al.  (9),  Yoshida  et  al.  (68),  Nishikawa  et  al.  (44,  46), 
and  Strickland  et  al  (64)  as  shown  in  Table  5. 

A  direct  comparison  of  FP  rates  is  made  even  more  difficult  due  to  the  lack  of 
standardized  images  areas,  and  image  area  has  shown  to  be  an  important  factor  in  FP 
detection  rates  (18).  As  shown  in  Table  5,  several  researchers  use  images  half  the  area  of 
those  used  in  this  research.  Additionally,  as  shown  in  Table  3,  if  subimages  are  used,  this 
research  could  claim  FP  rates  as  low  as  0.2295  if  image  area  is  ignored.  Nevertheless,  the 
detection  scheme’s  FP  rate  of  5.75  per  full-breast  image  is  comparable  with  reported  rates  of 
1  to  5  FPs  per  (sub)image  (9,  44,  46,  64,  68)  as  shown  in  Table  5.  (Because  Strickland  et  al 
used  50  ^m/pixel  images,  this  corresponds  to  a  10x10  cm^  region  on  a  film  mammogram, 
and  the  region  of  film  covered  in  this  research  corresponds  to  a  20x10  cm^  area.)  Moreover, 
if  image  area  is  used  to  adjust  FP  rates  for  comparison  purposes,  then  the  Ochoa  et  al 
detection  scheme’s  FP  rate  is  the  lowest  of  the  results  shown  in  Table  5. 
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Additionally,  although  other  investigators  do  not  provide  the  following  level  of  detail 
when  reporting  detection  results,  the  following  malignant  and  benign  detection  rates  were 
identified: 

•  Malignant  detection  rates 

—  94.1%  (16/17)  TP  case  detection  rate 

—  81.6%  (31/38)  TP  ROI  detection  rate 

—  80.0%  (28/35)  TP  full-breast  image  detection  rate 

—  5.6  FPs  per  full-breast  image  (20x10  cm^,  12  bit,  100  fira) 

•  Benign  detection  rates 

—  100%  (11/11)  TP  case  detection  rate 

—  92.0%  (23/25)  TP  ROI  detection  rate 

—  95.2%  (20/21)  TP  full-breast  image  detection  rate 

—  6.0  FPs  per  full-breast  image  (20x10  cm^,  12  bit,  100  //m) 

Because  other  investigators  do  not  report  this  level  of  detail  when  describing  their  detection 
schemes,  a  direct  comparison  of  these  detection  rates  is  not  an  option.  However,  the 
fairly  large  number  of  images  used  to  estimate  malignant  or  benign  TP  and  FP  rates  is 
comparable  with  the  size  of  the  entire  database  used  by  Yoshida  et  al.  (68),  Nishikawa 
et  al.  (44,  46),  and  Strickland  et  al  (64)  as  shown  in  Table  5. 

4.4  Summary 

The  database  of  56  digitized  (12  bit,  100  nm)  full-breast  (20x10  cm^)  film  mammo¬ 
grams  contained  63  biopsy-truthed  ROIs  over  28  cases.  This  technique  demonstrated  a  TP 
case  detection  rate  of  96.43  percent  (27/28)  and  TP  ROI  (54/63)  and  TP  image  (48/56) 
detection  rates  of  85.7  percent  with  5.75  FPs  per  full-breast  image. 
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V.  Conclusion 


5. 1  Summary 

This  research  proposes  the  following  methodology  for  clustered  microcalcification 
detection.  First,  preprocess  the  digitized  film  mammogram  to  reduce  digitization  noise. 
Second,  spatially  filter  the  image  with  a  difference  of  Gaussiajis  (DoG)  kernel.  To  detect 
potential  microcalcifications,  segment  the  filtered  image  using  global  and  local  threshold¬ 
ing.  Next,  cluster  and  index  these  detections  into  regions  of  interest  (ROIs).  Identify  ROIs 
on  the  digitized  image  (or  hardcopy  printout)  for  final  radiologic  diagnosis.  Finally,  to 
improve  detection  rates,  globally  optimize  detection  parameters  using  a  genetic  algorithm 
(GA),  then  locally  optimize  using  the  simplex  method. 

As  mentioned  earlier,  a  direct  comparison  of  detection  results  against  other  CADx 
detection  schemes  (9,  44,  46,  64,  68)  cannot  be  done  since  CADx  researchers  use  different 
databases.  However,  as  shown  in  Table  5,  the  85.7  percent  ROI  TP  detection  rate  found 
using  this  detection  scheme  is  comparable  with  reported  ROI  TP  detection  rates  of  84  to 
93  percent  (9,  44,  46,  64,  68).  Although  the  performance  of  a  CAD  scheme  is  difficult 
to  reliably  estimate  (45),  the  number  of  images  used  to  determine  TP  and  FP  rates  was 
larger  than  those  used  by  Chan  et  al.  (9),  Yoshida  et  al.  (68),  Nishikawa  et  al.  (44,  46), 
and  Strickland  et  al  (64)  as  shown  in  Table  5. 

A  direct  comparison  of  FP  rates  is  made  even  more  diflBcult  due  to  the  lack  of 
standardized  images  areas,  and  image  area  has  shown  to  be  an  important  factor  in  FP 
detection  rates  (18).  As  shown  in  Table  5,  several  researchers  use  images  half  the  area  of 
those  used  in  this  research.  Additionally,  as  shown  in  Table  3,  if  subimages  are  used,  this 
research  could  claim  FP  rates  as  low  as  0.2295  if  image  area  is  ignored.  Nevertheless,  the 
detection  scheme’s  FP  rate  of  5.75  per  full-breast  image  is  comparable  with  reported  rates 
of  1  to  5  FPs  per  (sub)image  (9,  44,  46,  64,  68)  as  shown  in  Table  5.  Moreover,  if  image 
area  is  used  to  adjust  FP  rates  for  comparison  purposes,  then  the  Ochoa  et  al  detection 
scheme’s  FP  rate  is  the  lowest  of  the  results  shown  in  Table  5. 
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5.2  Contributions 


The  following  list  summarizes  significant  contributions  made  during  the  course  of 
this  research: 

•  A  complete,  simple  to  use,  end-to-end  microcalcification  detection  system  was  devel¬ 
oped  which  demonstrated  a  96.4%  case  TP  rate,  and  an  85.7%  TP  rate  with  a  5.75 
FPs  rate  over  56  biopsy-truthed  full-breast  images  (20x10  cm^,  12  bit,  100  fj,m). 

•  This  system  is  the  first  reported  clustered  micro  calcification  detection  technique  de¬ 
veloped  using  optimization  methods  rather  than  mostly  ad-hoc  modification  of  de¬ 
tection  parameters. 

•  If  the  average  targetsize  is  known,  the  DoG  methodology  presented  here  may  easily 
be  used  in  other  target  recognition  schemes  because  the  algorithm  is  independent  of 
image  area  and  has  been  shown  to  be  insensitive  to  large  variations  in  image  quality; 
i.e.,  subtle  targets  in  noisy  backgrounds. 

•  The  relatively  straightforward  optimization  approach  presented  here  may  also  be  used 
in  other  target  recognition  research  to  efficiently  determine  optimimum  detection 
parameter  settings. 

•  Three  of  AFIT’s  most  recent  CADx  schemes  (10,  38,  48)  were  integrated  into  the 
first  AFIT  CADx  system  for  batch-detecting  radiographic  ROIs  in  digitized  film 
mammograms. 

•  During  the  course  of  this  work,  an  image  histogram  equalization  utility  was  de¬ 
veloped  and  included  in  Abrahamson’s  (1)  Pulse-Coupled  Neural  Network  (PCNN) 
magnetic  resonance  imaging  (MRI)  research.  The  equalization  step  was  included  in 
a  preprocessing  stage  to  allow  a  PCNN  to  accurately  segment  brain  MRIs. 

5.3  Conclusions  and  Recommendations 

The  objective  of  this  thesis  was  to  design  an  automated  microcalcification  detection 
system  which  may  be  used  by  radiologists  as  an  aid  in  mammogram  interpretation.  This 
objective  was  met  as  indicicated  by  a  TP  case  detection  rate  of  96.43  percent  (27/28)  and 
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TP  ROI  (54/63)  and  TP  image  (48/56)  detection  rates  of  85.7  percent  with  5.75  FPs  per 
full-breast  image  (20x10  cm^,  12  bit,  100  /xm). 

To  evaluate  the  clinical  impact  of  this  system,  a  clinical  study  should  be  performed 
on  this  research,  as  well  as  the  other  CADx  techniques  (10,  38,  48)  included  in  the  foremen- 
tioned  batch-capable  system.  Also,  to  reduce  the  FP  rate,  future  AFIT  CADx  researchers 
should  investigate  the  utility  of  fusing  the  results  of  this  research  with  other  clustered 
microcalcification  detection  methods.  Alternatively,  augmenting  this  system  with  a  neu¬ 
ral  network  using  features  derived  from  detected  ROIs  may  also  help  reduce  the  FP  rate. 
Finally,  the  optimized  DoG  filtering  approach  developed  during  the  course  of  this  research 
has  strong  potential  for  detecting  small  targets  in  other  cluttered  imagery,  such  as  detecting 
SCUD  missile  launchers  in  “noisy”  long-range  infrared  images. 
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Appendix  A.  Database  Information 


A.l  Overview 

In  screen-film  mammography,  the  film  is  used  for  detection,  storage,  and  display. 
(15).  Since  full-area  detectors  suitable  for  digital  mammography  are  still  in  the  experi¬ 
mental  stage,  an  alternative  way  of  acquiring  mammograms  in  digital  form  is  to  digitize 
a  film  mammogram  (67).  This  research  was  performed  using  digitized  film  mammograms 
and  subimages  annotated  by  a  physician  trained  in  radiology. 

Specifically,  the  database  of  56  digitized  (12  bit,  100  pm)  full-breast  (20x10  cm^) 
film  mammograms  contained  63  biopsy-truthed  clustered  microcalcification  ROIs  over  28 
cases.  All  mammograms  included  a  pathology  report  indicating  location  and  diagnosis  of 
biopsied  regions.  The  data  set  included  approximately  2  films  per  case  with  the  following 
distribution:  1  case  with  4  films,  2  cases  with  only  1  film,  and  the  remaining  25  cases 
with  2  films.  Full-breast  image  optimization  would  have  been  extremely  costly;  approx¬ 
imately  30  to  60  minutes  per  digitized  film  on  a  Sun  Ultra  or  Sun  Sparc20.  Thus,  52 
256x256  and  9  512x512  subimages  were  extracted  from  these  images  to  facilitate  algorith¬ 
mic  development.  The  61  subimages  contained  all  63  biopsy-truthed  ROIs,  and  a  physician 
annotated  abnormalities  identified  in  corresponding  pathology  reports.  These  subimages 
considerably  reduced  algorithmic  development  times,  and  reduced  computer  memory  re¬ 
quirements.  Additional  information  on  digitization  and  database  management  is  located 
in  the  Appendix. 

The  following  sections  were  written  by  Capt.  William  E.  Polakowski,  Capt.  Edward 
M.  Ochoa,  and  Jill  Leighner. 
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A.S  Digitization  of  film  mammograms 


Film  Mammogram  Digitization  Instructions 

I.  Power  up  the  Big  Mac  (Mac  lifx)  in  the  X-ray  filing  room,  (if  necessary) 

II.  To  calibrate  the  digitizer. 

A.  Note  that  the  digitizer  only  needs  to  be  calibrated  once  per  digitization  session. 

B.  From  the  large  folder  in  the  center  of  the  screen  with  three  icons  on  it,  click  on  Film 
Digitizer. 

C.  Note  that  to  enter  data  into  Mac  text  entry  windows,  you  must  hit  the  TAB 
key  instead  of  return. 

D.  Login  as  A  FIT  then  press  TAB. 

E.  Type  in  the  password  then  press  RETURN. 

F.  From  the  button  bar  on  top  of  the  screen,  select  “Special  then  Calibrate”. 

G.  Obtain  any  6, 14”  x  1 7”  films  from  underneath  the  digitizer. 

H.  Load  the  6  films  into  the  digitizer,  which  is  located  to  the  right  of  the  Big  Mac. 

I.  To  begin  the  calibration  procedure,  select  “OK”. 

III.  Obtain  and  use  the  1 6  hour  bar  code.  (Do  this  while  calibrating.) 

A.  Login  to  the  small  Mac  (Mac  llsi)  located  to  the  left  of  the  Big  Mac.  There  should  be  a  large 
folder  in  the  center  of  the  screen  with  two  icons  on  it. 

B  Select  Paris 

C.  Below  the  box  stating  Patients  Name,  type  “AFIT,  TEST”.  Then  select  Find  Patient. 

D.  From  the  top  of  the  Paris  window  click  on  icon  called  “Exams  for  patient”.  (An  exam 
is  one  digitization  session.) 

E.  Under  the  list  of  exams  for  AFIT,  TEST  select  any  exam.  (It  does  not 
matter  which  exam  is  selected  during  one  digitization  session.) 

F.  On  the  top  of  the  screen,  on  the  button  bar,  select  “Exams  then  Duplicate”. 

1 .  (If  necessary),  in  the  center  right  portion  of  the  screen,  in  the  Requesting  Ward/Clinic 
box,  type  “MAMMOGRAPHY”. 

G.  At  the  bottom  of  the  screen,  click  the  “New  Exam”  button. 

H.  Remove  the  barcode  slip  from  the  machine  to  the  left  of  the  small  Mac. 

I.  Select  “File  then  Quif.  Remember  the  barcode  slip  is  good  for  1 6  hours. 

IV.  To  digitize  the  films. 

A.  Scan  the  barcode  and  listen  for  the  beep  from  the  Big  Mac. 

B.  Prepare  mammography  films  you  intend  to  digitize. 

1 .  Make  sure  that  the  films  are  as  clean  as  allowable.  (If  OK,  wipe  off  any  grease  pencil  or 
smudge  marks,  and  if  possible  remove  stickers  -  they  tend  to  make  the  digitizer  refuse 
to  scan  the  film.) 

2.  Orient  the  films  as  the  radiologist  views  them. 

3.  Order  the  films  by  stacking  them  (sticker  side  down)  in  the  following  bottom  to  top 
order :  RCC,  LOG,  RMLO,  LMLO  (if  all  are  available). 

4.  Flip  the  stack  (top  toward  you)  so  the  RCC  view  is  now  on  the  top  of  the  stack  and 
rotate  the  top  of  the  films  away  from  you. 

5.  Remember,  the  digitizer  pulls  the  films  from  the  top  of  the  stack  to  the  bottom,  so  the 
final  order  of  digitization  will  be  RCC,  LCC,  RMLO,  LMLO. 

C.  Place  and  align  the  films  in  the  digitizer,  put  the  long  side  of  the  films  flush  against  the  right 
side  of  the  digitizer’s  auto  feeder. 

D.  To  begin  the  digitization  process,  on  the  Big  Mac,  in  the  Film  Digitizer  window,  click 
“Digitize”  then  click  “OK”. 

E.  After  the  films  have  completed  digitization,  click  “OK’. 

F.  To  end  digitization,  go  to  the  top  of  the  screen,  and  from  the  button  bar,  select  “File  then 
Quit”  from  the  pull  down  menu,  then  select  “Quit”  when  the  Login  window  reappears. 

G.  Return/Store  films  that  have  been  digitized. 


47 


A. 3  Digitized  film  storage  instructions 


Digitized  Mammogram  Storage  Instructions 

I.  To  connect  portable  hard  drive  for  digitized  film  downloading. 

A.  Follow  directions  carefully. 

B.  Locate  the  drive  adapter  cable  for  the  portable  Mac  drive. 

C.  To  connect  the  portable  Mac  drive  to  the  Big  Mac. 

1 .  From  the  main  window,  select  “Special  then  Shutdown”.  THE  SYSTEM 
MUST  BE  SHUT  DOWN  BEFORE  CONNECTING  THE  PORTABLE 
HARD  DRIVE. 

2.  To  attach  the  drive. 

a.  Remove  the  digitizer’s  black  connector  in  back  of  the  Big  Mac. 

b.  Attach  the  portable  drive’s  power  cord.  Make  sure  the  portable  drive’s 
power  switch  is  OFF.  (The  switch  in  the  back  of  the  drive  should  be  to 
the  empty  circle  position.) 

c.  Attach  the  portable  drive  to  the  Big  Mac’s  socket  the  digitizer  was  connected 
to.  Check  to  make  sure  that  the  keyboard  was  not  disconnected  during 
this  process. 

3.  To  activate  the  drive. 

a.  Note  that  the  next  few  steps  may  need  to  be  repeated  until  the  Big  Mac 
recognizes  the  portable  drive. 

b.  Turn  on  the  drive  (switch  In  the  back  of  the  drive).  The  light  in  front  will 
turn  off. 

c.  Press  the  Big  Mac  startup  key.  If  the  Big  Mac  fails  to  turn  on  check  that 
the  keyboard  was  not  disconnected  while  attaching  the  drive  or  push 

the  reset  switch  in  the  rear,  bottom  right  corner  of  the  back  of  the  Big  Mac. 

d.  As  the  Big  Mac  powers  up,  the  drive  will  turn  on,  then  off. 

e.  From  the  large  folder  in  the  middle  of  the  screen  select  “File  then 
Go  to  Finder”  and  enter  the  password. 

f.  To  use  each  disk. 

1 .  Insert  disk  (Be  sure  to  keep  track  of  the  ones  you  have  already  used.) 

2.  Open  the  numbered  icon  (click  on  It). 

3.  Clear  each  disk  so  it  can  hold  the  digitized  mammograms.  Select 
from  the  top  menu  “Special  then  Erase  Disk  then  Initialize” 

II.  To  get  image  information  and  download  images  to  disks. 

A.  From  the  large  folder  in  the  main  window,  double  click  on  FWD060  and  double 
click  on  Lite  Box. 

B.  Login  as  A  FIT,  then  press  TAB. 

C.  Type  the  password,  then  press  RETURN. 

D.  To  obtain  the  image  name  and  size  information. 

1 .  Go  to  the  Patient  Name  box  and  type  AFIT,  TEST  (no  return  key  is  required 
here,  It  will  automatically  find  the  AFIT  list  for  you). 

2.  Select  your  exam  (the  last  exam  at  the  bottom  of  the  list).  It  should  be  an  exam 
with  an  “I”  In  the  ONLINE  column.  The  “i”  means  that  the  process  is  still  active. 

3.  To  display  all  of  the  exam  images,  on  top  of  the  “All  exams”  window ,  click  on 
“Images  then  Display”. 

4.  Select  OK  when  prompted. 

5.  Click  on  the  first  image. 

6.  To  obtain  specific  image  information. 

a.  From  the  top  button  bar  select  “Image  then  Image  About”. 

b.  Record  filenames  (only  the  last  8  characters)  and  image  sizes,  the  size 
is  in  the  middle  of  the  list  and  the  filename  is  near  the  bottom. 

c.  Repeat  for  all  remaining  images. 

E.  To  download  images  to  portable  disks. 

1 .  Close  AFIT,  TEST  window. 

2.  Select  “File  then  Open”  (the  same  exam  should  already  be  selected). 

3.  Select  “Images  then  Verify”  to  change  “I”  to  “Y”. 


48 


4.  Select  “Images  then  Store  Exam”. 

5.  Make  sure  the  numbered  disk  is  already  selected. 

6.  Select  “Store”.  When  the  disk  is  full,  remove  the  disk  and  insert  a  new  one  (make 
sure  that  you  are  using  clear  disks).  Note  that  each  disk  holds  approximately  12-1 5 
mammograms. 

F.  To  exit  Lite  Box  select  “File  then  Quit”  from  the  top  button  bar,  DO  NOT  SAVE  CHANGES. 

III.  To  reset  the  Big  Mac  and  digitizer. 

A.  To  remove  the  remaining  150MB  disks  from  the  drive,  select  “Special  then  Eject  Disk”. 

B.  From  the  main  window  select  “Special  then  Shutdown”.  THE  SYSTEM  MUST  BE  SHUT 
DOWN  BEFORE  RECONNECTING  THE  DIGITIZER. 

C.  Carefully  reconnect  the  digitizer.  (Reverse  drive  connect  steps.) 

D.  Press  the  restart  button  and  then  the  startup  key. 

IV.  Store  the  drive  adapter  cable. 
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A, 4  transfer  instructions 


Digitized  Mammogram  Image  Transfer  Instructions 

I.  To  transfer  digitized  mammogram  files  to  AFIT. 

A.  Go  to  room  201 1 . 

B.  Locate  the  adapter  cable. 

C.  Shut  down  the  host  system. 

D.  Connect  the  portable  hard  drive. 

E.  Insert  disk. 

F.  Open  the  numbered  drive  window. 

G.  Select  local  exams  and  then  open  the  numbered  folder. 

H.  Now  rename  the  files,  removing  all  but  the  last  8  characters. 

I.  Open  Vulcan  and  select  Applications  folder. 

J.  Open  TCP/Connect  II  Folder  and  select  TCP/Connect  II. 

K.  Select  “Services”  on  the  button  bar,  then  “FTP”. 

L.  Click  on  “more  options”  and  unselect  “Use  Web-Style  FTP  client” 

M.  Set  the  host  name  to  pinna  and  set  username. 

N.  Click  Image  data  type. 

O.  Under  “Options”,  select  “Binary”  and  unselect  “Prompt  for  filenames”. 

P.  Set  the  directory  for  storage  (for  example  /home/pi nnal/bdata/Uncut). 

Q.  Click  on  Desktop,  select  the  files  and  click  on  “Copy”. 
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Breast  Cancer  Group  Filenaming  Guide 


Breakdown  of  a003a05m.dog:  a  003  a  0  5  m  .dog 

a 

Source 

a:  Wright-Patterson  AFB  Medical  Center,  WPAFB,  OH 
g:  Grandview  Hospital,  Dayton,  OH 
s:  Southview  Hospital,  Dayton,  OH 

f:  Franciscan  Health  Center,  Dayton  Campus,  Dayton,  OH 

003 

Case  Number 

Each  case  typically  contains  1-4  digitized  mammograms,  and 
corresponds  to  a  single  date 

a 

Orientation 

a:  RCC 
b:  LCC 
c:  RMLO 
d:  LMLO 

e:  RML 

f:  LML 

0 

Date  Number 

Dates  are  represented  in  increasing  order,  where  the  earliest 
is  represented  by  0  and  the  latest  is  the  highest  integer  value. 

5 

Subimage  Index 

If  not  0,  then  this  number  is  used  to  index  subimages  ex¬ 
tracted  from  each  digitized  mammogram. 

m 

Diagnosis 

m:  Malignant,  Cancer 
b:  Benign 

u:  Unknown,  no  biopsy 
n:  Not  available  yet 

.dog 

File  Type 

dog  -  Difference  of  Gaussians  filtered,  grayscale,  BIN 

hit  -  Hit  and  miss  filtered,  grayscale,  BIN 

mam  -  Original  mammogram,  grayscale,  BIN 

nul  -  Region  without  lesions,  256x256  pixels,  grayscale,  BIN 

pen  -  PCNN  processed  SUBIMAGE 

roi  -  SUBIMAGE,  256x256  pixels,  grayscale,  BIN 

ro2  -  SUBIMAGE,  512x512  pixels,  grayscale,  BIN 

tru  -  Radiologist  and/or  pathology  verified  truth,  ASCII 

txt  -  Comments  (corresponding  date,  etc.),  ASCII 

unc  -  Uncut  (full  digitized  mammogram),  grayscale,  BIN 

wav  -  Wavelet  processed,  grayscale,  BIN 

51 


Appendix  B.  Matlab  and  C  Code 

The  file  eucdet.mis  the  parent  Matlab  function  used  to  automatically  detect  clustered 
micro  calcifications.  The  CMEX  file,  findcluster.c  was  written  by  Capt.  Amy  Magnus  and 
modified  by  Capt.  Edward  Ochoa  to  compile  on  Solaris. 


GetFileNames.m 


function  [FullFileNames, JustFileNames]  =  GetFileNames (ListOf Files, DirOfFiles) ; 
y.  GETFILENAMES  Used  to  load  ascii  filenames  into  matlab. 

'/,  function  FullFileNames  =  GetFileNames  (ListOf  Files , DirOfFiles)  ; 
y.  ASSUMES:  All  filenames  have  the  same  number  of  characters 
y.  EX:  FullFileNames  =  GetFileNames(*filenames.txt  Whome/DataBaseDirO 
y,  EX:  FullFileNames  =  GetFileNames( ^filename s. txt *  ,0)  */,  if  in  cwd 

I 

%  by:  Capt.  Edward  M.  Ochoa,  GE0“96D 

FilelD  =  fopen(ListOfFiles , ^r* ) ; 
if  FilelD  ==  -1, 

error(sprintf( TERROR:  '/s  did  not  open  correctly  (GetFullFileNames .m) * ,ListOfFiles)) ; 
end; 

FileList=fread(FileID) ; 

FileLi8t=setstr(FileList O ; 

for  f ilenum=l :sum(isspace(FileLi8t)) 
if  isstr(DirOfFiles) 

if  "isdir(DirOfFiles) , 

err  or  (sprintf(»  ERROR:  '/.s  is  not  a  directory  (GetFullFileNames  .m)  », DirOfFiles)  )  ; 
end; 

[FileN,FileList]=8trtok(FileLi8t); 

FullFileNames (filenum,  :)=sprintf(’y.8/y.s’  , DirOf Files  ,FileN)  ; 

JustFileNamesCfilenum,  :)=sprintf  (’y,s>  ,FileN)  ; 

else 

[FullFileNames (filenum, : ) ,FileList]=strtok(FileLi8t) ; 


end 

FileNamelD  =  fopenCFullFileNames (filenum, :), *r’ ) ; 
if  FileNamelD  ==  -1, 

error (sprintf  (TERROR:  y,s  not  found\nCheck  DirOfFiles  (GetFullFileNames.m) \FullFileNames (filenum, :))) 
end; 

f close (FileNamelD) ; 

end 

fclose(FilelD) ; 


GetMammo .m 


function  mammo=GetMammo(fullfilename,NUMROWS,DISPLAYTF) ; 
y,  GETMAMMO  Used  to  load  a  breast  cancer  database  image, 
y,  mammo=GetMammo(fullfilename,NUMROWS,DISPLAYTF) ; 
y,  LOADS  FULL  MAMMO,  REGARDLESS  OF  SIZE 
y.  by:  Capt  Amy  Magnus,  Capt.  Edward  M.  Ochoa,  GE0-96D 

FileNotFound  =  0; 

fid  =  fopen(fullfilename,  *r*); 

if  fid  ==  -1 
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FileNotFound  =  1; 
end; 

*/  IF  FILE  IS  FOUND,  LOAD  DIGITIZED  MAMMOGRAM 
if  "FileNotFoTind 

[mammojcnt]  =  freadCf id,inf ,  *ushortO  ; 
numcols  =  cnt/NUMROWS; 

maitimo  =  r eshape(maimno,NUMROWS, numcols )  ; 

clear  cnt 
f closeCf id) ; 

else 

errorCsprintf  ( **/,8  not  found’ ,fullfilename))  ; 

end; 

if  DISPLAYTF 

mamodi sp=mammo (1:4: NUMROWS ,1:4: numcols ) ; 
colormap (gray (4095 ) ) 
image (mamodisp) 
axis  image 

title(8printf  ( *y.s’  ,fullf  ilename))  ; 
axis  off 

end 


PutMammo.m 


function  PutMammo(M,ful If ilename) 

%  PUTMAMMO  Used  to  write  out  an  image  using  the  breast  cancer  database  format. 
'/,  PutMammo(M,fullf ilename) 

'/  This  m-file  writes  the  matrix  M  to 

y,  the  file  ’  fullf ilename  ’  .  The  onus 

y,  is  on  the  user  to  remeber  M’s  dimensions. 

y.  by:  Capt  Amy  Magnus,  Capt.  Edward  M.  Ochoa,  GE0-96D 

fid=f open (fullf ilename,  ’w+’); 
fwrite(f id,M, ’ushort ’ ) ; 
f close(f id) ; 


ebuildpt .m 

function  detect  =  ebuildpt (xon, yon, numrows, numcols) ; 
y.  EBUILDPT  Given  a  set  of  pixel  coordinates,  construct  an  mask  image, 
y,  detect  =  ebuildpt  (xon, yon, numrows, numcols )  ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

detect=zeros (numrows, numcols) ; 

for  i=l: length (xon) 

detect (xon(i) ,yon(i) )=1 ; 

end 


ecenters .m 


function  [centers] =ecenters (inlist ,map) ; 
y*  ECENTERS  Compute  centroids  of  clusters, 
y,  [centers]=ecenters( inlist, map) ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0~96D 

centers=zeros (size (map)) ; 
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for  i=l : size (map, 1) ; 

xy=inlist (map(i , 1) :map(i ,2) , : ) ; 
if  size(xy,l)==l 
centers (i, : )=xy ; 
else 

centers  (i , :  )=roTind(mean(xy) )  ; 

end 

end 


echeckpts .m 

function  [go odpt s]® echeckpts (or igimg,pt coords) ; 
y,  ECHECKPTS  Used  to  verify  detections  were  made  within  breast 
'/.  [goodpts]=echeckpts(origimg,pt coords) ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

if  ptcoords(l,l)"=-l 

mask=esegmam(origimg) ; 
n\impts=s ize(pt coords  ,1) ; 
keep=ones (numpts , 1) ; 
for  i=l: numpts 

if  " mask (pt coords (i,l) ,ptcoords(i,2)) 
keep(i)=0; 
end 

end 

goodpts®pt coords (keep , : ) ; 
if  isempty(goodpt8) 
goodpts=[-l  -1] j 

end 

else 

goodpts=pt coords ; 

end 


ecrnois .m 

function  [outimg]=ecrnois(img) ; 

y,  ECRNOIS  Cross-shaped  hybrid  median  filter  an  image. 

%  [outimg]=ecmois(img) ; 

y*  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

winsize=3; 

[numrows ,numcols]®size(img) ; 
bufr  ®  (winsizo-l)/2; 

y.  MY  "MIRROR"  PADDING  TO  ELIMINATE  EDGE  EFFECTS  IN  CODE 
img=emirrpad(img,winsize) ; 
y,  INITIALIZE  OUTPUT  IMAGE; 
outimg  =  zeros  (nujnrows,numcols) ; 


y.  ELIMINATE  NOISE 

grpv=zeros(6,l) ; 

temp®0; 

j=0; 

i=0; 

minval=0 ; 
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minind=0 ; 


for  m=l:niimrows 
mb=m+bufr ; 
for  n=l:niimcols 
nb  =  n+bufr; 


y,  GET  +  PIXELS 
grpv(l)=img(mb,nb-l) ; 
grpv(2)=img(mb,nb) ; 
grpv(3)=img(mb,nb+l) ; 
grpv(4)=img(mb-l,nb) ; 
grpv(5)=img(mb+l,nb) ; 

for  i=l:3 

ininval=grpv(i)  j 
minind^i ; 
for  j=i+l:6 

if  grpv(j)<miiival 
ininval=grpv  ( j ) ; 
miiiind=j ; 

end 

end 

if  (minval<grpv(i)) 
teitip=grpv(i)  ; 
grpv(i)=minval; 
grpv  (itiinind)  =temp ; 
end 

end 

outiing(m,n)*grpv(3) ; 

end 


end 


edetclst .m 


function  [inli at , map] =odet cist (list, radius ,numpts) ; 
y,  EDETCLST  Used  to  detect  clusters  of  pts 
y,  [inlist ,  map]  =edot  cist  (list,  radius,  numpts) ; 
y»  by:  Capt.  Edvard  M.  Ochoa,  GE0-96D 

pairs  =  epairs (list , radius ) ; 
inlist=C]  ; 
map=  []  ; 

if  "is empty (pairs) 

[ind,indmap]=ere coir (pairs) ; 
startr=l; 

for  i=l :sizo(indmap,l) 

inds-ind(indmap(i,l) :indmap(i,2)) ; 
clstr=list (inds , : ) ; 
numcpts=size(clstr , 1) ; 
if  numcpts>=numpts 

inlist= [inlist ;clstr] ; 
endr=startr+numcpts“l ; 
map=[map;startr  endr] ; 
startr=endr+l; 

end 
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end 


end 


if  is empty (inlist) 
inlist=[“l  -1] ; 
map^inlist ; 
end 


edogk.m 


function  [dogkarnel,Gl,G2]=edogk (target size) ; 

'/  EDOGK  Used  to  construct  the  DoG  kernel  given  targetsize. 
%  [dogkernel,Gl,G2]=edogk(targetsize) ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

if  (targetsize>33)  |  (targetsize<=0) 

disp(* Sorry,  limited  to  0  <  targetsize  <  33  pixels’); 
return 
end 

sig  =  [.391  .625  1  1.6  2.56  4.096  6.536  10.486  16.777]; 
Nsize  =  [3  5  9  11  15  21  29  43  55] ; 

s2  =  targetsize/2; 
sl=s2/1.6; 

y.  Find  optimum  N  (kernel  size) 

i=l; 

while  sig(i)<s2 
i=i+l; 

end 

N=Nsize(i) ; 

G1  =  f speciaK’gaussian’ ,N,sl)  ; 

G2  =  f speciaK’gaussian’ ,N, 82)  ; 
dogkernel  =  G1  -  G2; 


edrawROI .m 


function  edrawROI (inlist , map ,numrows ,numcols ,drawPtsTF) 
y.  EDRAWROI  Used  to  highlight  ROIs 
y,  edrawROI  (inlist  ,map,numrows  ,numcols  ,drawPtsTF) 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

for  i*l:size(map,l) 

index=map(i,l) :map(i,2) ; 

M= [inlist (index, 2)  inlist (index, 1)] ; 
if  i==l 

if  det(cov(M)) 

ellipse (mean (M) ,inv(cov(M)) ,3.5) ; 
else 

axis(’ij  ’) ; 

plot(M(:,l),M(:,2),’r+’); 

end 

axis(’ij ’) ; 
hold  on; 

if  drawPtsTF  &  (det(cov(M))) 
plot(M(:,l),M(:,2),’b+’); 
end 

axi8([l  numcols  1  numrows]); 
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axis image*) ; 
else 

if  drawPtsTF  ft  (det (cov(M))) 
plot(M(:.l),M(:.2),*b+*); 
end 

if  det(cov(M)) 

ellipse (mean (M) ,inv(cov(M)) ,3.5); 
else 

axis(*ij*); 

plot(M(:,l),M(:,2),*r+*); 

end 

end 

end 


efilt .m 


function  [outimg]  =  ef i It (img,numbits, kernel ) ; 

y,  EFILT  Handy  utility  for  convolving  (filtering)  an  image  with  a  kernel. 
y»  [outimg]  =  efilt (img,numbits , kernel)  ; 

'/,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

[numrow8,numcols]=size(img) ; 

N=size (kernel,!) ; 
bufr  =  (N~l)/2; 

y,  MY  "MIRROR"  PADDING  TO  ELIMINATE  EDGE  EFFECTS  IN  CODE 
img=emirrpad(img,N) ; 

y.  RESCALE  THE  IMAGE  TO  INTENSITY  VALUES  [0,1]; 
img=erescale(img) ; 

y,  CONVOLVE  THE  KERNEL  WITH  THE  IMAGE 
out img=conv2(img, kernel, ’valid*) ; 

y,  RESET  THE  OUTPUT  IMAGE  TO  THE  COLORDEPTH  DESIGNATED  BY  NUMB  ITS 
outimg=round(  (2 '‘numbits - 1)  *erescale (outimg) )  ; 


eget truth. m 

function  t rut h=egettruth ( SOURCED IR,IMGNAME, class) ; 
y.  EGETTRUTH  Utility  for  loading  pathology  truth, 
y,  truths5egettruth(S0URCEDIR,IMGNAME, class) ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

TRUTHFULLFILENAME=[SOURCEDIR  IMGNAME  ’ .tru’] ; 
evaKsprintf  (’load  ‘/.s;  ’  ,TRUTHFULLFILENAME) ) 
evaKsprintf  ( ’truth=y,s;  ’ , IMGNAME)) ; 
targetindox=f ind( truth ( : ,5)==class) ; 
if  "is empty (tar get index) 

truth=t ruth (target index, : ) ; 
else 

truth=-l; 

end 


eglth2.m 


function  glthr=eglth2(f img,NUMBITS,globp,segwin8ize,k,ktop) ; 

y,  EGLTH  Implements  global  and  local  thresholding  for  segmentation/detection, 

y,  glthr=eglth2(fimg,NUMBITS,globp,segwinsize,k,ktop) ; 
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y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 


gl=ehistgr(fimg,NUMBITS,globp) ; 

glthr=elthresh2 (f img ,NUMBITS ,segwinsize ,k ,gl ,ktop) ; 
glthr=bwmorph(glthr, ^bridge O ;  %  BRIDGE  1  PIXEL“WIDE  GAPS 


ehist .m 


function  hgram  =  ehist (img,nnmbits) ; 

'/.  EHIST  Utility  for  computing  an  image  histogram. 

y,  hgram  ®  ehist  (img,numbits) ; 

y,  by:  Capt.  Edwsird  M.  Ochoa,  GE0~96D 

ind=0; 

img=round ( img ( : ) ) ; 
limg=size(img,l) ; 

hgram=zeros(  realpow(2,numbits)  ,  1  ); 
for  i=l:limg 
ind=img(i)+l; 
hgram (ind)=hgr am (ind)+l ; 
end 


ehistgr .m 

function  [gr , imhs t ] =ehi st gr ( img , NUMB ITS , percentile) 

y,  EHISTGR  Used  to  segment  an  image  using  a  percentile  of  an  image  histogram, 
y,  gr=ehi  stgr  (  img ,  NUMB  ITS ,  percent  ile) 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

imhs t=ehist (img, NUMBITS) ; 
numcolr8=length(imhst) ; 
chst=cumsum(imhst) ; 
cdf=chst/chst(numcolrs) ; 
gr=min(f indCcdf >=percentile)) ; 


ellipse. m 


function  ellipse(m,Q ,k) ; 

'/  ELLIPSE  Used  to  plot  an  ellipse, 
y,  ellipse (m,Q,k) ; 

X 

y,  Plots  the  ellipse  described  by 

y. 

y,  {x:  (x-m)’Q(x-m)  =  rsq} 

y. 

y,  where 

y,  m  =  mean  of  data  points,  (vector) 

y,  Q  =  inverse  of  covariance  of  data  points .  (matrix) 

y,  k  ®  number  of  std  devs  to  include  in  ellipse.  (3  works  well) 

y. 

y*  USAGE:  ellipse(mean([xpts  ypts]  )  ,  inv(cov(  [xpts  ypt8])),k) 
y,  by:  Maj  Ruck 

m  =  m( : ) ; 
pts  *  2B6; 

circle^cmplx  =  exp(j*2*pi/(pts-l)*[0:pts-l] ) ; 
circle  =  [real(circle«cmplx) ;  imag(circle_cmplx)] ; 
ellip  =  sqrt(k)*Q"(-0.6)*circle  +  m*ones (1 ,pts) ; 


ecolorl=*k’ ; 
lwidthl=2 ; 
stylel=*- * ; 
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ecolor2=*r* ; 
lwidth2=l ; 
style2=»-»; 

hl*plot(ellip(l, :) ,ellip(2, :)) ; 
set (hi , * LineWidth ’ , Iwidthl ) ; 
set (hi , *  Color  * , ecolorl) ; 
set(hl,  *LineStyle\styl6l) 
hold  on 

h2=plot(ellip(l, :) ,ellip(2, :)) ; 
set(h2, * LineWidth* ,lwidth2) ; 
set (h2 , *  Color  * , ecolor2) ; 
set (h2 , *  LineStyle  * , style2) 


elthresh2.m 


function  [outimg]=elthresh2(img,numbits,winsize,LTH,MINTHRESH,LTHT0P) ; 

ELTHRESH  Implements  local  and  global  thresholding 
y*  [outimg]  =elthresh2(img,numbits ,vinsize ,LTH,MINTHRESH,LTHTOP)  ; 

*/»  by:  Capt .  Edward  M.  Ochoa>  GE0~96D 

MAXGRAY=2"numbits-l; 
y.  BEGIN  THRESHOLD  TEST 

[numrows ,numcols] =size(img) ; 
bufr  =  (winsize-l)/2; 

y,  MY  "MIRROR”  PADDING  TO  ELIMINATE  EDGE  EFFECTS  IN  CODE 
img=emirrpad(img,winsize) ; 


y.  INITIALIZE  OUTPUT  IMAGE; 
outimg  =  zeros (numrows, numco Is ) ; 
y.  COMPUTE  FIRST  ROW  WINDOW  *S  TOTALS 


tot*0; 

01dtots=0 ; 
tot0fsquares=0; 

OldtotOf squares=0 ; 
winsize2=realpow(winsize,2) ; 

for  m=l:winsize 
for  n=l:winsize 

tot  =  tot  +  img(m,n); 

tot  Of  squares  =  tot  Of  squares  +  img(m,n)''2; 

end 

end 

y.  INITIALIZE  FIRST  ROW  WINDOW’S  RESULTS 
01dtot=tot ; 

OldtotOf squares  =  totOf squares ; 

y,  COMPUTE  LOCAL  MEAN  AND  STANDARD  DEVIATION 
mn  =  tot/winsize2; 

sdev  =  realpow(  (totOfsquares  ~  tot“2/winsize2)/(winsize2-l) ,  1/2  ); 

y,  COMPUTE  FIRST  THRESHOLD  TEST 
8eed=img(bufr+l ,bufr+l) ; 

if  (seed>(mn  +  LTH*sdev))  ft  (80ed<(mn  +  LTHTOP*sdev))  ft  (seed>  MINTHRESH) 
outimgd ,  1)"1; 
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end 


y,  COMPUTE  THRESHOLD  TEST  FOR  REST  OF  IMG 

for  m=l:miinrows 
inb=in+buf  r ; 
tot=01dtot ; 

totOf  squeiressOldtotOfsquares ; 
for  n=2:miincols 
nb  =  n+bufr; 
for  indox=-bufr :bufr 

oldpix  *  img(mb+index,nb“bufr-l) ; 
newpix  =  img (mb+ index ,nb+bufr) ; 
tot  =  tot  -  oldpix  +  newpix; 

totOf squares  =  totOf squares  “  oldpix'‘2  +  newpix'‘2; 

end 

y.  FOR  COLUMNS  2:NUMC0LS,  COMPUTE  THRESHOLD  TEST 
mn  =  tot/winsize2; 

sdev  =  realpow(  (totOf  squares  “  tot'*2/winsize2)/(winsize2-l) ,  1/2  ); 
seed=img(mb,nb) ; 

if  (seedXmn  +  LTH*sdev))  ft  (soed<(nm  +  LTHTOP*sdev))  ft  (seed>  MINTHRESH) 
outimg(m,n)*l; 

end 


end 

y,  SINCE  AT  LAST  COL,  THEN  FOR  NEXT  ROW,  COLUMN  1, 
y,  COMPUTE  THRESHOLD  TEST 

if  m'snumrows 
nb  =  bufr  +  1; 
tot=01dtot; 

totOf squares=01dtot0f squares ; 
for  index=-bufr:bufr 

oldpix  ®  img(mb-bufr,nb+index) ; 
newpix  =  img (mb+bufr+l,nb+ index ) ; 
tot  =  tot  -  oldpix  +  newpix; 

totOf  squares  =  totOf  squares  -  oldpix'‘2  +  newpix'‘2; 

end 

Oldtot  =  tot; 

Oldtot Of squares  =  totOf squares; 
mn  =  tot/winsize2; 

sdev  s=  realpowC  (totOf  squares  -  tot''2/winsize2)/(winsize2-l) ,  1/2  ); 
seed=img(mb+l,bufr+l) ; 

if  (seed>(mn  +  LTH*sdev))  ft  (seed<(mn  +  LTHTOP’i'sdev))  ft  (seed>  MINTHRESH) 
out img (m+1 , 1) =1 ; 

end 


end 

end 


emailmsg.m 

function  emailmsg(emailaddress , message) ; 
y,  EMAILMSG  Used  to  send  email  message  to  user, 
y,  emailmsg(emailaddress,  mess  age)  ; 
y,  by:  Capt.  Edward  M,  Ochoa,  GE0-96D 
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UCMDl=sprintf  (’echo  >  /tmp/emailmsg .txt ’  .message) ; 

UCMD2=sprintf  (’echo  ""s  ’/.s"  »  /tmp/emai lmsg.txt’ .message)  ; 
UCMD3=sprintf (’echo  »  /tmp/emailmsg.txt ’ ) ; 

UCMD4=sprintf(’mail  */,s  <  /tmp/emailmsg.txt  ’  .emailaddress) ; 
UCMD5=sprintf (’ !rm  /tmp/emailmsg.txt’); 

UCMDS=[’!  ’  UCMDl  ’;’  UCMD2  ’;’  UCMD3  ’;’  UCMD4 

eval(UCMDS) 

eval(UCMD6) 


emakeodd.m 


function  oddnum  =  emakeodd(nimi) ; 

y,  EMAKEODD  Used  to  convert  input  to  closest  odd  integer  (ex:  60  ->  61) 
*/,  oddnum  =  emakeodd(num) ; 

'/«  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

rnum=round(num) ; 
ndif f=num-rnum ; 

isoddTF=rem(rnum,2) ; 

if  "isoddTF 

orignum=rnum+num ; 
upnum=rnum+l ; 
downnum=rnum~l ; 

if  abs (upnum-orignum) <=ab8 (downnum-orignum) 
oddnum=upnum ; 
else 

oddnum=downnum ; 

end 

else 

oddnum=rnum; 

end 


emirrpad . m 


function  img=emirrpad(img,winsize) 
y,  EMIRRPAD  Used  to  mirror  pad  image, 
y,  img= emirrpad (img, wins ize) 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

[numrows ,numcols]=size(img) ; 
bufr  =  (winsize-l)/2; 


y,  MY  ’’MIRROR"  PADDING  TO  ELIMINATE  EDGE  EFFECTS  IN  CODE 

tbuf  =  flipud(img(2:l+bufr, :)) ; 

bbuf  =  f 1 ipud ( img (numrows “bufr: numrows - 1, :)) ; 

corner  =  ones(bufr .bufr) ; 

Ibuf  »  [corner ; f liplr (img ( : , 2 : 1+buf r) ) ; corner]  ; 

rbuf  ®  [corner ;f liplr (img( : ,numcols-bufr :numcols-l)) ; corner] ; 

img  =  [Ibuf  [tbuf ; img; bbuf]  rbuf]; 

tic  =  f lipud(f liplr (img (1+buf r:winsize, 1+buf r : winsize)) ) ; 
trc  =  flipud(f liplr (img (1+buf r:winsize, numcols :numcols+bufr))) ; 
blc  =  flipud(f liplr (img (numrows: numrows+bufr, 1+buf r: winsize))) ; 
brc  =  flipud(f liplr(img(numrows :numrows+bufr, numcols :numcols+bufr))) ; 
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img ( 1 : 1+buf r , 1 : 1+buf r ) =tlc ; 

img ( 1 : l+buf r , numcols+buf r : niuncols+winsize- 1) =tr c ; 
img (numrows+buf r : numrows+winsize-1 , 1 : 1+buf r) =blc ; 

img(numrows+bufr : numrows+winsize-1, numcols+buf r:numcols+winsize-l)=brc ; 


epairs . m 

function  pairs  =  epairs (list , radius ) ; 
*/,  pairs  =  epairs  (list , radius) ; 

'/,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

numpts=size(list,l) ; 
pairs=C] ; 
for  i=l:numpts 
ptl*list(i, :) ; 
for  j=i+l:numpts 
pt2=list(j , ; 
if  dist(ptl,pt2)<=radius 
pairs  =  [pairs;!  j] ; 
end 

end 

end 


erecolr .m 


function  [list , map] =erecolr (crash) ; 

55  ERECOLR  Used  to  identify  the  number  of  unique  colors  given  crash  list 
55  [list  ,map]  «erecolr( crash) ; 

55  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

list=[]  ; 
map=  □  ; 

blist=l; 

tmp=crash; 
pairs=[]  ; 
i=0; 

keepchockingTF=l ; 
pairs=tmp(l , : ) ; 
while  keepcheckingTF 
for  j=l: size (pairs, 2) 
rows=[]  ; 

if  "isempty (tmp) 

[rows , cols] =find(pairs(j)==tmp) ; 

end 

if  "isempty(rows) 

rows=eunique(rows) ; 
newp=tmp (rows , : ) ; 
newp^eunique (newp) ; 
pairs=eunique(  [pairs  newp] ) ; 
tmp(rows  ,:)  =  □; 
end 

end 

if  isempty(tmp) 

keepcheckingTF=0 ; 

list=[list;  pairs(:)]; 
elist=blist+size(pairs ( : )  ,  D-l ; 
map=[map;bliBt  elist] ; 
blist=elist+l; 

elseif  ''any(evsamo(pairs,tmp( :))) 
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list=[list;  pairsC:)]; 
elist=blist+sizo(pair8 ( : ) , 1)-1 ; 
maps [map jblist  elist]  ; 
blist=elist+l; 

pairs=tinp(l , : ) ; 
end 

end 


erescale .m 


function  imgout  =  ere scale ( img) ; 

*/.  ERESCALE  Used  to  rescale  any  image  to  the  [0,1]  range. 

*/,  imgout  =  erescale  (img) ; 

y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

lo  =  min(min(img)) ; 
hi  =  max(max(img)) ; 

imgout=(img-lo)/(hi-lo) ; 


erocem.m 


function  [truepos,  falsepos,  goodroi] =erocem (truth, inlist , map) ; 
y,  EROCEM  Utility  for  identifying  whether  a  CLUSTER  was  detected, 
y  [truepos,  falsepos,  goodroi]=erocem(truth, inlist , map) ; 
y*  by:  Capt.  Edward  M.  Ochoa,  GEO-OGD 

y,  A  cluster  is  considered  a  TP  if  any  of  the  following  conditions  exist 

y,  1:  ALL  cluster  pixels  were  within  a  truth  box. 

y,  2:  At  least  MINUETS  were  found  within  a  corresponding 

y,  truth  box  and  the  centroid  of  the  cluster  is  within 

y,  the  truth  box  perimeter. 

y,  3:  At  least  MINUETS  were  found  within  a  corresponding 

y.  truth  box  and  the  centroid  of  the  cluster  is  within  MINUIST 

y,  pixels  from  the  truth  box  perimeter. 

y. 

y,  Note:  (1)“(3)  may  be  combined  into  the  following  statement: 

y.  At  least  MINUETS  were  found  within  a  corresponding  truth  box 

y.  and  the  cluster's  centroid  is  within  MINUIST  pixels  from  the  truth  box 

MINUIST=50;  */,  6mm  for  lOOum  digitized  images 
MINUETS=1; 

if  truth(l,l)"=-l 

truepos=0; 

good=ones(8ize(map,l) ,1) ; 
goodroi=ones(size(map,l) ,1) ; 

centroids«ecenters (inlist ,map) ; 
for  i=l; size (map, 1) 
good(i)=0; 

index=map(i ,1) :map(i ,2) ; 
centroidsscentroids  (i , : )  * ; 

y.  FOR  EACH  CLUSTER,  COUNT  THE  NUMBER  OF  UETS  WITHIN  A  TRUTH  BOX 
for  k=index 

for  j=l:size(truth,l) 
rlo=truth(j ,2) ; 
rhi=truth(j ,4) ; 
row=inlist(k,l) ; 
rin=(rlo<=row)ft(row<=rhi) ; 
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clo=truth(j ,1) ; 
chi=truth(j ,3) ; 
col=inlist(k,2) ; 
cin=(clo<=col)&(col<=chi) ; 
inboxTF  =  (rin  &  cin) ; 
if  inboxTF 

good(i)=good(i)+l; 

end 

end 

end 

numgood=good(i) ; 
good(i) =good(i) >MINDETS-1 ; 

%  IF  MINDETS  WITHIN  A  TRUTH  BOX  IS  SATISFIED  AND  ALL  DETS  ARE  NOT 
'/,  WITHIN  A  TRUTH  BOX,  CHECK  IF  CENTROID  WITHIN  TRUTH  BOX 
if  good(i)  &  (niimgood"=size(index,2)) 
crow=centroid(l) ; 
ccol=centroid(2) ; 
for  j=l: size (truth, 1) 
rlo=truth(j ,2) ; 
rhi=truth(j ,4) ; 
rin2=(rlo<=crou)ft(crow<=rhi) ; 
clo=truth(j ,1) ; 
chi=truth(j ,3) ; 
cin2*(clo<=ccol)&(ccol<*=chi) ; 
centroidinboxTF  =  (rin2  &  cin2) ; 

end 

if  centroidinboxTF 
good(i)=l; 

y,  IF  MINDETS  WITHIN  A  TRUTH  BOX  IS  SATISFIED  AND  ALL  DETS  ARE  NOT 
y,  WITHIN  A  TRUTH  BOX,  AND  CENTROID  NOT  IN  BOX,  CHECK  MINDIST 
y,  OF  CLUSTER  CENTROID 
else 

good(i)=0; 

for  j*l : size (truth, 1) 
rlo=truth(j ,2) ; 
rhi=truth( j ,4) ; 
clo=truth(j ,1) ; 
chi=truth(j ,3) ; 

boxpts  «  etrubox(rlo,rhi,clo,chi) ;  */,  PERIMETER  PIXELS 

centroidv=centroid*ones (1 , size (boxpts , 1) ) ; 

y,  dists®dist  (boxpts  ,controidv) ; 

dists=edist (boxpts , centroidv) ; 

inindistTF=2my  (any  (di8ts<MINDIST)  )  ; 

if  mindistTF 

good(i)=good(i)  |  1; 
end 

end 

end 

end 

if  ''good(i) 

goodroi(i,l)=0; 

end 

end 

truepos=sum(goodroi==l) >0 ; 
f alsopos«sum(goodroi==0) ; 

else 

truepos=0; 

falsepo8=size(map,l) ; 
goodroi=-l*ones(size (map,l)) ; 
end 
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esegmam.m 

function  mainmask  -  esegmamCmam) ; 

%  ESEGMAM  Used  to  segment  breast  mass  from  full  mammogram. 

mammask  =  esegmam(mam) ; 

'/  by:  Capt .  Edward  M.  Ochoa >  GE0-96D 

[numrows  ,numcols3  -sizeCmeim) ; 

•/,  SUBSAMPLE  TO  SIMPLIFY  SEGMENTATION 
mams  =mam (1:16: numrows ,1:16: numcols ) ; 

'/.  HYBRID  MEDIAN  FILTER  TWICE  TO  MERGE  REGIONS 
mams =ecrnois (mams) ; 
mams =ecrno is (mams) ; 

I  COMPUTE  AND  RESCALE  VARIANCE  IMAGE  (TO  ENHANCE  BREAST  EDGE) 
mamv=erescale (log (1+evarop (mams ,12,9))); 

'/,  ADD  VARIANCE  IMAGE  TO  ORIGINAL  (SIMPLIFIES  SEGMENTATION) 
mamsv=erescale(mamv+ere scale (mams) ) ; 

t  USE  FLOYD-STEINBERG  DITHERING  TO  CONVERT  TO  BW 
mamma8k“dither(mam8v) ; 

*/,  APPLY  MORPHOLOGICAL  OPERATORS  TO  SEGMENT  BREAST 
mammas k=bwmorph (mammask, » majority* ,4) ; 
mammask=bwmorph (mammask , *open* ,2) ; 
mammask=bwmorph (mammask , *  erode  * , 3) ; 

y,  REMOVE  LABEL  IF  NECESSARY 
manmiask=findcluster (mammask) ; 
maxlabel=max (mammask ( : ) ) ; 
maxc=l ; 

for  i=l:maxlabel-l 

if  sum  (mainmask  ( :  )==i)<sum  (mammask  ( :  )==(i+l) ) 
maxc=i+l; 
end 

end 

mammask=mammask==maxc ; 

[numr,numc]*size(mEuns) ; 

y,  ORIENT  IF  NECESSARY  (TO  SIMPLIFY  FILL) 
lhs=sum(mams ( : , 1) ) ; 
rhs=8um(mams ( : ,numc) ) ; 
if  lhs>rhs 
ImamTF— 1 ; 
else 

lmamTF=0 ; 

end 

if  "ImamTF 

mammask=f lip Ir (mammask) ; 
end 

y,  HORIZONTAL  FILL 
col8um=sum(mammask) ; 
ind-min(f ind(max(colsum)==colsum)) ; 
fillcol=mammask(: ,ind); 
for  i=l:ind 

mammask( : ,i)=mamma8k(: ,i)  |  fillcol; 
end 
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I  VERTICAL  FILL 
rows\im=siim(mainma8k* ) ; 
begrow=min(find(rowsum)) ; 
endrow=max(f  indCrowsum) ) ; 

topr=min(f  ind(rowsum*=max(rowsimi(f  ind(rowsum''=0)  )  )  )  )  ; 
botr=max(f  ind(rowsuiri==Tiiax(rowsum(f  indCrowsum^sO))  ) ) ) ; 

sumr*max(rowsuiti) ; 
for  i=begrow : topr 

if  (rowsumCiXsumr)  &  (rowsiim(i)''=0) 
siimr=rowsiiin(i) ; 
toprind=i ; 

end 

end 

toprf ill=mamma8k(toprind, : ) ; 
for  i=begrow ; topr 

mammas k ( i,  :)=mainmask(i, :)  I  toprfill; 
end 

siimr*max(rowsiim) ; 
for  i=botr :endrow 

if  (rowsum(i)<sumr)  ft  (rowsTim(i)''=0) 

8iimr=rowsum(i) ; 
botrind=i ; 
end 

end 

botrf  ill=mainmask(botrind, : )  ; 
for  i=botr : endrow 

mammaskCi ,  :)=mainmask(i , :)  I  botrf  ill; 
end 

'/,  CONVERT  BACK  TO  FULL  SIZE  MASK 
if  "ImamTF 

mammask=fliplr(mammask) ; 
end 

mammask=imr e 8 ize( mammas k, [numrows  numcols] , ^nearest  O ; 
mammask=roimd(mammask) ; 


esepname.m 

function  [name, ext, f name, nrow8]«e8epname(FULLNAMES,isnaun0TF) ; 

*/.  ESEPNAME  Handy  string  manipulator  used  to  break  down  filenames,  and  ID  num  rows 
*/»  [n8Uiie,ext  ,f name, nr ows]=esepnamo (FULLN AMES, isnameTF) ; 

*/.  by:  Capt,  Edward  M.  Ochoa,  GE0-96D 

f  name=  []  ; 
name=  []  ; 
ext~[]  ; 

if  "isnameTF 

revfullname=fliplr(FULLNAMES) ; 

for  i=l:8ize(FULLNAMES,l) 

[n,rmdor]=strtok(revfullname(i, :),*/*); 
f name= [f name ; n] ; 

end 

fname«f liplr(fname) ; 
else 

fname=FULLNAMES; 

end 

for  i*l : size (f name, 1) 
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[n,rmder]=strtok(fname(i, . 0  ; 
name = [name ;n] ; 

ext  = [ext ; rmder (2 : s ize (rmder , 2 ) ) ] ; 

end 

for  i=l:8ize(ext,l) 

if  3s:=sTiin(ext(i,  :)==*mam* ) 
nrows(i,l)®2048; 
elseif  3==siim(ext(i,  :)®=*roi’ ) 
nrows(i,l)=256; 

elseif  3==8um(ext(i, :)==’ro2’) 
nrows(i,l)®612; 
end 

end 


eshrink.m 


f\inction  [pts  ,pt coords]  =  eshrink(detmask)  ; 

y.  ESHRINK  Used  to  reduce  detections  to  single  points.  Facilitates  clustering. 
'/t  [pts  ,pt coords]  =  eshrink(detmask)  ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 


if  sum(8um(detmask) )==0 

pts=zeros(size(detmask) ) ; 
else 

[numrows ,numcols]=8ize(detma8k) ; 
detmask=f indclusterCdetmask) ; 

y,  Get  histogram  on  ’values*  image 
h=ehist(detmask,8) ; 
obj^colors=find(h''=0) ; 
nc^length (obj . colors )- 1 ; 
ptcoord8=one8(nc,2) ; 

if  (nc"=0) 
for  i=l:nc 

[x ,y] =f ind(detmask®®i) ; 
if  size(x,l)®®l 
xy= [x  y] ; 
else 

xy=round(mean(  [x  y]  )  )  ; 

end 

ptcoordsCi, : )=xy; 

end 

end 

pts®ebuildpt (pt coords ( : ,1) ,ptcoords(: ,2) , numrows ,numco Is) ; 


etrubox.m 


function  boxpts  =  etrubox(rlo,rhi,clo,chi) ; 

y,  ETRUBOX  Utility  for  converting  truth  indicies  to  perimeter  pts 

y,  boxpts  =  etrubox (truth) ; 

y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

topcol®  [do :  chi]  ; 
toprow=rlo*ones(size(topcol)) ; 

rgtrow=rlo+l :rhi ; 
rgtcol=chi*ones(size(rgtrow)) ; 
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botcol=chi-l:-l :clo; 
botrow='rhi*ones(siz6(botcol))  ; 


If trow= [rhi“l : -1 : rlo+1] ; 
lftcol=clo*ones (size (If trow)) ; 

boxpts-Ctoprow^  topcol* ; . . . 
rgtrow*  rgtcol* ; . .  ► 
botrow^  botcol^ ; . . . 

If trow ^  Iftcol*]; 


eucdet .m 


function  . . . 

[inlist , map, nimirows,numcols,roiTF,TPROI,FPROI,truth]=. . . 

eucdet (  PARAMS , . . . 

IMGNAMEEXT, . . , 

DoStatsTF, . . . 

DisplayTF) ; 

•/,  EUCDET  Detect  Clustered  Microcalcifications:  Capt.  Edward  M.  Ochoa,  GE0-96D 

y. 

y,  [inlist, map, nujnrows,numcols,roiTF,TPROI,FPROI,truth]=. . . 
y,  eucdet (  PARAMS,...  <—  DETECTION  PARAMETERS  (described  below) 
y.  IMGNAMEEXT,... 

y,  DoStatsTF,...  <~-  Compute  ROC  stats?  (yes=l,  no=0) 
y,  DisplayTF);  <—  Display  processing  results?  (yes=l,  no=0) 

y. 

y,  DESCRIPTION  OF  DETECTION  PARAMS  =  [pi  p2  p3  p4  p5  p6  p7]  AND  DEFAULT  SETTINGS 
y.  pi  -  DOGT,  targetsize  (1.1995) 

'/  p2  -  K,  multiple  of  local  std  for  local  thresholding  (4.1226) 

y,  p3  -  GLOBP,  global  histogram  percentile  cutoff  (0.9991) 

y*  p4  ■“  DETWINSZ,  odd  window  size  used  in  local  thresholding  (69) 

y,  p5  -  CLSTRRAD,  max  dist  between  neighbor  microcalcs  (61.3019) 

y.  p6  •*  NUMCALCS,  min  niimber  of  calcs  to  trigger  a  cluster  detection  (3) 

y,  p7  “  KTOP,  max  multiple  of  local  std  for  local  thresholding  (7) 

y 

y,  EXAMPLE  ~  HOW  TO  USE  DEFAULT  SETTINGS: 

y*  [inlist  ,map,numrows,numcols]=eucdet(n  ,  *a036a01m.roi*  ,1»1)  J 
y,  [inlist  ,map,numrows,numcols]=eucdet(n  ,  * a036a00m.mam*  ,1,1); 

% 

I  EXAMPLE  “  HOW  TO  DISPLAY  DETECTION  RESULTS: 
y,  edr  awRO  I  (inlist,  map,  numrows,numcols,l) ; 

y. 

y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

y. - - ALGORITHM  SETTINGS . . - 

y,  DEFAULT  PARAMETER  SETTINGS 
if  isempty (PARAMS) 

PARAMS= [1.1996  4.1226  0.9991  69  61.3019  3  7]; 

end 

y.  IMG  INFO 

[IMGNAME , EXT , NAME , IMGNUMROWS] =esepname (IMGNAMEEXT , 0) ; 

ProcMamTF* (sum((l*EXT®sones (size (EXT,  1)  ,l)*»mamO  0==3)  * ; 

y.  DIR,  EXT,  AND  IMG  DEPTH  INFO 
if  ProcMamTF 
SOURCEDIR®*  * ; 

MAMTRUDIR® ’ /home/hawkey e3/96d/aochoa/Thes is /MC.code/MAMTRU/ ’ ; 
else 

SOURCEDIR=Vhome/pinnal/bdata/wpafbh/ROIs/'  ; 
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end 


NUMBITS=12; 

t  BLOCK  PROCESSING  SIZE  (REDUCES  MEMORY  REQUIREMENTS  AND  SPEEDS  UP  PROCESSING) 
BLKSZ=266; 

BLKSZR=BLKSZ; 

BLKSZC=BLKSZ; 

%  DoG  TARGETSIZE  and  filter  KERNEL 
D0GT=PARAMS(1); 

[KERNEL, dogkl,dogk2]  =  edogk(DOGT); 

y.  DETECTION  (SEGMENTATION)  PARAMETERS 
K=PARAMS(2) ; 

KT0P=PARAMS(7); 

GL0BP=PARAMS(3) ; 

DETWINSZ=emakeodd(PARAMS(4)) ; 

y,  CLUSTERING  PARAMETERS 
CLSTRRAD=PARAMS(5); 

NUMCALCS=roimd(PARAMS(6)) ; 

t  -  ALGORITHM  BEGINS  - 

y,  LOAD  IMG 

img=GetMammo( [SOURCED IR  IMGNAME  EXT  ] , IMGNUMROWS , 0) ; 

y,  img=GetMainmo(IMGNAMEEXT, IMGNUMROWS, 0);  */.  if  IMGNAMEEXT  always  includes  fullpath 
[niimrows  ,mimcols]=size(img) ; 

if  ProcMamTF  &  (siiin(img( :  ,numcols))>siiin(iitig(  :  ,niimcols)) ) 

FlipTF=l ; 
img=fliplr(img) ; 
else 

FlipTF=0; 

end 

if  DisplayTF  I  ProcMamTF 
origiing=img; 
else 

origimg=-l; 

end 

if  DisplayTF 

NUMGRAY=2-*NUMBITS; 

figure 

imsliow((NUMGRAY-l)*erescale(img) ,gray(NUMGRAY)) 
title (sprintf(  ^y,s^  ,NAME)) 
orient  tall 

end 

X . — - - - - - - - — . - . 

y.  CLEAN  NOISE 

img=blkproc(img, [BLKSZR  BLKSZC] , * ecrnois ’ ) ; 

if  DisplayTF 
figure 

imshow((NUMGRAY-l)*erescale(img) ,gray(NUMGRAY) ) 
t it le  ( sprint f  (» y,8 .  ncl » ,  IMGNAME)  ) 
orient  tall 

end 

X - - - - - - — 

X  DoG  FILTER 
P1=NUMBITS; 

P2=KERNEL; 

img=blkproc(img , [BLKSZR  BLKSZC] , » ef ilt  * ,P1 ,P2) ; 
if  DisplayTF 
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figure 

ks iz6=size( KERNEL, 1) ; 

subplot (2,1,1) ,m6sh(KERNEL) ,axis  square 

axis  ij 

axis  (Cl  ksize  1  ksize  min(  [dogkK : )  *  dogk2(:)0)  max(  [dogkl(  :)  ^  dogk2(:)0)]) 

title(’DoG  filter  coefficients’) 

subplot (2,1, 2) ,freqz2 (KERNEL, 32,32) ; ,  axis  square 

title (’DoG  filter  Frequency  Response’) 

textsc ( . 6 , . 97 , sprint f  ( ’target size=*/,l . 2f  ’  ,D0GT) ) 

orient  tall 

figure 

iitishow  (  img ,  gray  (NUMGRAY  )  ) 
title(sprintf ( ’’/.s.dog’  ,IMGNAME)) 
orient  tall 

end 

I . . 

'/  USE  GLOBAL/LOCAL  THRESHOLDING  TO  SEGMENT/DETECT  POTENTIAL  TARGETS 
P1=NUMBITS; 

P2=GL0BP; 

P3«DETWINSZ; 

P4=K; 

P5=KT0P; 

img=blkproc(img, [BLKSZR  BLKSZC] , ’eglth2’ ,P1 ,P2 ,P3 ,P4,P5) ; 

if  DisplayTF 
figure 

imagesc(img) ,  colormap( 1-gray (2)) , axis  image 
t it le  ( spr intf  ( ’  '/.s .  thr » ,  IMGNAME) ) 
orient  tall 

end 

I - - - 

y,  SHRINK  DETECTIONS  TO  SINGLE  PIXELS/POINTS 
img=blkproc( img, [BLKSZR  BLKSZC] , ’eshrink’ ) ; 

if  DisplayTF 
figure 

images c( img) ,  colormap( 1-gray (2)) , axis  image 
title (sprintf(  ’Jls.pts’  , IMGNAME)) 
orient  tall 

end 

t - - - - - - 

*/,  ELIMINATE  FULL  MAMMO  EDGE  DETECTIONS  (BASED  ON  KERNEL  SIZE) 
if  ProcMamTF 

N=size (KERNEL,!); 
bufr=(N-l)/2; 

img(l:bufr, :)=zeros(bufr,numcol8) ; 
img(numrows-bufr+l rnumrows , :)=zeros(bufr,numcols) ; 
img( : , 1 :bufr)=zero8(numrows ,bufr) ; 
img( : ,numcol8-bufr+l :numcols)=z6ro8 (numrows ,bufr) ; 

end 

if  FlipTF 

img=fliplr(img) ; 

end 


I - 

y,  ID  DETECTION/POINT  LOCATIONS 
[xon,yon] =f ind(img==l) ; 

if  isempty(xon) , 
ptcoord8*[-l  -1]; 
else 

ptcoords=[xon  yon]  ; 
end 
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% . 

I  CHECK  IF  POINTS  WITHIN  BREAST 
if  ProcMamTF 

[pt coords] =Gcheckpts ( or igimg,pt coords) ; 

end 

I - - - — - - - -  - - - 

'/.  USE  CLUSTERING  CRITERIA  TO  DETECT  ROIS  AND  INDEX  PTS  WITHIN  ROIS 
[inlist ,map] =Gdet cist (pt coords ,CLSTRRAD,NUMCALCS) ; 


if  DoStatsTF 

%  COMPUTE  ROI  DETECTION  STATISTICS 

class“2; 

if  "ProcMamTF 

truth=GgettrTith(SOURCEDIR, IMGNAME , class)  ; 
else 

t ruth=egett mth (MAMTRUDIR , IMGNAME , clas  s ) ; 

end 

if  inlist(l,l)"=-l 
if  truth" =-l 

[TPROI,  FPROI,  roiTF]=Grocem(truth,inlist,map) ; 
else 

TPR0I=O; 

FPROI=size(map,l) ; 
roiTF=zeros(size(map,l) ,1) ; 
end 
else 

TPROI-0; 

FPR0I=0; 

roiTF=“ones(sizG(map,l) ,1) ; 

end 

else 

truth=-l; 

TPR0I=“1; 

FPR0I=-1; 

roiTF=-l; 

end 

%  - . - . . . . . - - - 

'/.  IF  DESIRED,  DISPLAY  RESULTS 
if  DisplayTF 
figure 

if  inlist (1, !)"=-! 

edrawRO I (inlist ,map,numrows ,numcol8 , 1) ; 
title (sprintf(*ROI(s)  in  '/s^NAME)) 
if  TPR0I"=-1 

xlabeKsprintf  (»TPR0I*y.d  FPROI^Kd* , TPROI, FPROI)) 

end 

else 

title (sprintfC* No  ROI(s)  in  '/.s*,NAME)) 
if  TPR0I"==-1 

xlabeKsprintf  (>TPR0I=y.d  FPR0I=y,d’ , TPROI, FPROI)) 

end 

end 

if  DoStatsTF 
axis  ij 

axis([l  numrows  1  numcols] ) 
showro i t rut h ( IMGNAME , 2 ) 

end 

orient  tall 
figure 

imshow((NUMGRAY-l)*erescale(origimg) ,gray(NUMGRAY)) 
hold  on 

if  inlistCl, !)"=-! 

edrawRO I ( inlist ,map , numrows , numcols , 0) ; 
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title (sprintf(^ROI(s)  in  */,s*,NAME)) 
if  TPRDI^=-1 

xlabeKsprintf  (>TPROI='/.d  FPROI=y,d'  ,TPROI,FPROI)) 
end 
else 

title (sprintfC* No  ROI(s)  in  */,s’,NAME)) 
if  TPR0I^=“1 

xlabeKsprintf  (*TPR0I=*/.d  FPR0I=*/d'  ,TPR0I,FPR0I)) 
end 

end 

if  DoStatsTF 

showroitruthClMGNAME ,2) 

end 

orient  tall 

end 

t - ALGORITHM  DONE - - — 


eunique . m 

function  uniquels  =  eunique(v); 

*/,  uniqels  =  eunique (v); 

y,  by:  Capt,  Edward  M.  Ochoa,  GE0“96D 

tmp=v( :) ; 

[nrows ,ncols]=8ize(v) ; 

uniquels  =[]  ; 
i=0; 

keepcheckingTF=l ; 
while  keepcheckingTF 
i=i+l ; 

ind=f  ind(v(i)==titip) ; 
totun=size(ind,l)  ; 
tinp(ind)  =  [] ; 

if  totun(l)  I  (totun==nrows) 
uniquels= [uniquels  v(i)] ; 
end 

if  isempty(tmp) 

keepcheckingTF=0 ; 

end 

end 

tmiquels=sort (uniquels) ; 


evaorop.iti 


function  [outiiiig]=evarop(img,numbits ,winsize) ; 

y,  EVAROP  Computes  the  standard  deviation  image  given  neighborhood  size  (odd) 
y,  [out img]=evELrop(img,numbits, wins ize)  ; 
y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 

MAXGRAY*roalpow(2,numbits)-l ; 

y.  BEGIN  VAR  ALGORITHM 

[numrows ,numcols] »size(img) ; 
bufr  =  (winsize”l)/2; 

y.  MY  "MIRROR”  PADDING  TO  ELIMINATE  EDGE  EFFECTS  IN  CODE 


img=emirrpad(img, wins ize) ; 
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t  INITIALIZE  OUTPUT  IMAGE; 


outimg  “  zeros (numrowsjiiiunco Is ) ; 

I  COMPUTE  FIRST  ROW  WINDOW »S  TOTALS 

tot=0; 

01dtot=0; 

tot0fsquares=0; 

0 Idt  ot  Of  squar es=0 ; 
wiiisize2=realpow(winsize,2) ; 

for  in=l:winsize 
for  ii=l:winsizG 
tot  =  tot  + 

totOf  squares  “  totOf  squares  +  realpow(img(in,n)  ,2)  ; 

end 

end 

*/.  INITIALIZE  FIRST  ROW  WINDOW  »S  RESULTS 
01dtot*tot ; 

OldtotOf squares  =  totOf squares ; 
y,  COMPUTE  LOCAL  STANDARD  DEVIATION 

sdev  =  realpowC  (totOf squares  -  tot*2/winsize2)/(winsize2-l) ,  1/2  ); 

y.  COMPUTE  FIRST  LOCAL  STANDARD  DEVIATION 
outimgCl ,l)=sdev; 

y.  COMPUTE  LOCAL  STANDARD  DEVIATION  FOR  REST  OF  IMG 

for  ms=l:numrows 
mb=m+bufr; 
tot=01dtot ; 

totOf squares=01dtot0f squares ; 
for  n=2:numcols 
nb  =  n+bufr; 
for  index=~bufr :bufr 

oldpix  =  img(mb+index,nb-bufr-i) ; 
newpix  =  img(mb+index,nb+bufr) ; 
tot  =  tot  -  oldpix  +  newpix; 

totOfsquares  =  totOf squares  -  realpow (oldpix, 2)  +  realpow (newpix, 2) ; 
end 

y,  FOR  COLUMNS  2:NUMC0LS,  COMPUTE  LOCAL  STANDARD  DEVIATION 

sdev  =  realpow(  (totOfsquares  -  tot“2/winsize2)/(winsize2-l) ,  1/2  ); 

outimg (m , n) =sdev ; 


end 

y.  SINCE  AT  LAST  COL,  THEN  FOR  NEXT  ROW,  COLUMN  1, 
y,  COMPUTE  LOCAL  STANDARD  DEVIATION 

if  m"=numrows 
nb  =  bufr  +  1; 
tot=01dtot ; 

totOf squares=01dtot Of squares ; 
for  index= -bufr : bufr 

oldpix  =  img(mb-bufr,nb+index) ; 
newpix  =  iing(mb+bufr+l,nb+index) ; 
tot  =  tot  -  oldpix  +  newpix; 
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totOf squares  =  totOlsquares  -  realpow(oldpix,2)  +  realpow(newpix,2) ; 

end 

Oldtot  =  tot; 

Oldtot Of squares  =  totOf squares; 

sdev  =  realpow(  (totOf squares  -  tot"2/winsize2)/(winsize2-l) ,  1/2  ); 
out  irag  (m+ 1,1)  =sde  v ; 

end 

end 

out img=M AXGRA  Y . * ere s  c  ale ( out img) ; 


evsame .m 


function  vs ame=evsaiiie(v, value ) ; 

y,  EVSAME  Used  to  identify  elements  common  to  both  vectors 
y,  vsame=6Vsamo (v, value) ; 
y,  by:  Capt.  Edward  M.  Ochoa 

vsame=zeros (sizo(v) ) ; 

for  i=l:length(v) 

for  j«l: length (value) 

vsame(i)=(v(i)==value(j))  I  vsame(i); 
end 

end 


myprint .m 


f unc t i on  myprint (f igs tart , f igend , har dcopyTF , cr e at ef i leTF ) 
y,  MYPRINT  A  handy  utility  to  print  figures 
y,  myprint  (figst art  ,f igend  ,hardcopyTF ,  createf  ileTF) 

y. 

y,  by:  Capt.  Edward  M.  Ochoa,  GE0-96D 
for  k=figstart:f igend 

if  (k<10) ,  f ignoeps=[’f igO’  int2str(k)  *.eps^]; 

else  fignoeps=[^f ig’  int2str(k)  ’.eps*];  end 
eval([’figure(*  int2str(k)  O*]) 
eval(C ’print  -deps  ’  fignoeps]) 
if  hardcopyTF,  eval( [’ Iqmslpr  ’  fignoeps]);  end 
if  "createf ileTF,  eval([Mrm  ’  fignoeps]);  end 

end 


showmam . m 


function  img= showmam (IMGNAMEEXT,displayTF) ; 
y,  SHOWMAM  Used  to  load/display  mammo  image 
y,  img=showmam(IMGNAMEEXT,di8playTF) ; 
y,  by:  Capt,  Edward  M.  Ochoa,  GE0-96D 

y.  IMG  INFO 

[IMGNAME , EXT , NAME , IMGNUMROWS] =esepname (IMGNAMEEXT , 1) ; 

Pr 0 cM8unTF= (sum( (l*EXT==ones (size (EXT, 1) ,1)*’ mam’) ’ )==3) ’ ; 

y.  DIR,  EXT,  AND  IMG  DEPTH  INFO 
if  ProcMamTF 

SOURCEDIR= ’ /home/pinnal/bdata/wpaf bh/ ’ ; 
else 
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SOURCED rR=  *  /home/piiinal/bdata/wpaf bh/ROIs/ » ; 

end 

NUMBITS=12; 

'/.  LOAD  IMG 

img=GetMainmo  ( [SOURCED  IR  IMGNAME  ».»  EXT  ]  ,IMGNUMR0WS,0) ; 

if  displayTF 
■/.  DISPLAY  IMG 

MAXGRAY=realpow (2 , NUMBITS ) - 1 ; 

imshow ( (MAXGRAY-1) ♦erescaleCimg) ,gray (MAXGRAY) ) 
title (sprintf(  ,NAME)) 
orient  tall 

end 


showroitruth.m 


function  showroitruth(ROINAME, CLASS) ; 

'/.  SHOWROITRUTH  Used  to  display  pathology  truth  box  given  class 
y,  showroitruthCROINAME, CLASS) ; 
y,  by  Capt.  Edward  M.  Ochoa,  GE0-96D 

y,SOURCEDIR=  *  /home/hawkoye3/96d/eochoa/Thesis/MC_code/MAMTRU/ '  ; 
SOURCED IR=’ /home/pinnal/bdata/wpafbh/ROIs/ ^ ; 
TRUTHFULLFILENAME=[SOURCEDIR  ROINAME  ».tru*]; 
eval(sprintf(*load  y.s;  MRUTHFULLFILENAME)) 
evaKsprintf  (*truth=y»s;  * , ROINAME))  ; 
showtruth (truth , CLASS ) 


showtruth.m 


function  not rut hTF= showtruth (truth , class ) 

y,  SHOWTRUTH  Used  to  display  pathology  truth  box  given  class  of  interest 
y,  no  t  rut  hTF=sho  wt  ruth  (truth ,  class) 
y,  by:  Capt.  Edweurd  M.  Ochoa,  GE0-96D 

if  truthd, !)''=-! 

hold  on 

targetindex^find (truth ( : ,6)== class) ; 

for  index=l : length(targetindex) ; 

target=truth(targetindex( index) ,1:4) ; 
xstart=target(2) ; 
ystart=target(i) ; 
xend  =tearget(4); 
yend  =target(3); 

xplot=[xstart  xstart  xend  xend  xstart]  -  .B; 
yplot«[ystart  yend  yend  ystart  ystart]  -  .5; 
hl=plot (yplot ,xplot , ) ; 
set(hl,*LineWidth>,2); 
if  index==l 
hold  on 
end 

plot (yplot , xplot , *  b  O 
end 

end 
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textsc.m 


function  H  =  tGxt8c(x,y ,txt) ; 

*/,  TEXTSC  places  text  in  screen  coordinates  and  places 
'/  a  title  at  the  top  of  the  figure. 

y. 

y.  H  =  TEXTSC (X,Y, TXT)  places  the  text  string,  TXT 
y,  at  the  noritialized  coordinates  X  and  Y .  His  the 
y,  handle  to  the  text  object. 

y. 

y,  H  =  TEXTSC  (TXT,  *  title*)  places  a  title  at  the  top 
y.  of  the  figure  window.  This  is  useful  when  you 
'/  want  to  place  a  single  title  above  multiple 
y,  subplots. 

y. 

y,  TEXTSC  creates  an  invisible  AXES  which  occupies 
y,  the  entire  FIGURE.  The  units  of  the  AXES  are 
y,  normalized  (range  from  0  to  1) .  TEXTSC  checks 
y,  all  the  children  of  the  current  FIGURE  to  determine 
y.  if  an  AXES  exist  which  meets  these  criteria  already 
y,  exist.  If  one  does,  then  it  places  the  text  relative 
y,  to  that  AXES . 

y. 

%  NOTE:  Requires  MATLAB  4.2  or  above. 


y,  Written  by  John  L.  Galenski  III 
y.  All  Rights  Reserved  January  21,  1994 
y.  LDM031695jlg 

y,  Basic  error  checking 
if  nargin  <  2 

error (* TEXTSC  requires  at  least  2  inputs*) 
end 

y,  Check  to  see  if  AXES  already  exist 
ch  =  get (gcf ,* Children’ ) ; 

ax  =  f indob j (gcf , ’Type’ , ’axes* , ’Tag* , ’TEXTSC*) ; 
if  isempty(ax) 

ax  =  axes ( ’Units ’, ’Normal’ , ’Position’ , [0  0  1  1],  ... 

’Visible’ .’Off ’, ’Tag’, ’TEXTSC’); 

else 

axes (ax) 
end 

y.  Place  the  text 

if  nargin  =-  2  &  isstr(x)  &  strcmp(lower(y) , ’title’ )  '/  Subplot  title 

txt  =  x; 

X  =  .6; 

tmp  *=  text  ( ’Units ’, ’normal’ , ’String’ , ’tmp’ , ’Position’ ,  [0  0  0]); 
ext  =  get (tmp , ’ Extent ’ ) ;  * 
delete (tmp) 

H  =  ext (4) ; 
y  =  1  -  .60*H; 
end 

h  =  text (x,y, txt, ’VerticalAlignment ’, ’Middle’ ,  ... 

’ HorizontalAlignment ’ , ’ Center ’ ) ; 

y,  Make  the  original  AXES  current 
if  ''isempty(ch) 

set (gcf , ’CurrentAxes’ ,ch(l)) 
end 

y,  Check  for  output 
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if  nargout  ==  1 
H  =  h; 
end 


f indcluster.c 

MEX  file  written  by  Capt.  Amy  Magnus. 

To  compile  on  Solaris,  modified  by  Capt.  Edward  M.  Ochoa 


#include  <stdio.h> 

#inclnde  <stdlib.h> 

#includo  <math.h> 

#inclnde  "mex.h” 

typedef  struct  { 
int  rows; 
int  cols; 
double  *ptr; 

}  matlab.MATRIX; 


static 

#ifdef  ..STDC_. 
void  color .cluster ( 
int  X, 
int  y, 

long  * color, 
matlab.MATRIX  Mask 
) 

#else 

color.cluster (x ,y , color , Mask) 
int  X,  y; 
long  *color; 
matlab.MATRIX  Mask; 

#endif 

{ 

int  i,j,  xi,yj; 

Mask.ptr [x*Mask .rows  +y]  =  (double)  *color; 

for  (i=“l;  i<2;  i++)-C 
xi  =  X  +  i; 

if  ((xi  >=  0)&&(xi  <  Mask.cols)){ 
for  (j=-l;  j<2;  j++){ 

yj  =  y  + 

if  ((yj  >=  0)&&(yj  <  Mask. rows) ){ 

if  (  Mask. ptr[xi*Ma8k. rows  +yj]  ==  -1  ){ 
color.cluster (xi,yj , color, Mask) ; 

} 


} 

} 

} 

} 

return ; 

} 

static 

#ifdef  __STDC.. 
void  find.clusters( 
matlab.MATRIX  Mask, 
matlab.MATRIX  IN, 
double  * threshold 
) 

#else 
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f ind.clusters (Mask , IN , threshold) 
matlab.MATRIX  IN , Mask ; 
double  ^threshold; 

#endif 

long  * color, zero; 
int  X,  y; 

for  (x=0;  x<Ma8k.col8;  x++){ 
for  (y=0;  y<Mask.rows;  y++){ 

Mask. ptrCx*Mask. rows  +y]  =  IN .ptr [x*IN.rows  +y]  >=  ^threshold  ?  -1  :  0; 

} 

} 

/♦initialize  colors/ 

zero  =  0; 
color  =  ftzero; 

/♦  find  the  first  pixel  of  cluster 
increment  color 

call  recursive  "color. cluster"  ♦/ 

for  (x=0;  x<Mask.cols;  x++){ 
for  (y=0;  y<Mask.rows;  y++)-C 

if  (Mask. ptr [x+Mask. rows  +y]  ==  -1){ 

♦color  +=  1; 

color.cluster(x,y , color , Mask) ; 

} 

} 

} 

return ; 

} 

#ifdef  __STDC__ 
void  mexFunction( 

int  nlhs , 

Matrix  ♦plhs [] , 
int  nrhs , 

Matrix  ♦prhs  [] 

) 

#else 

mexFunct ion (nlhs ,  plhs,  nrhs,  prhs) 
int  nlhs,  nrhs; 

Matrix  ♦plhs  []  ,  ♦prhs  []  ; 

#endif 

matlab.MATRIX  IN,  OUT; 
double  ♦threshold; 
double  default .threshold; 

if  ((nrhs  <  1)  | |  (nrhs  >  2)){ 

mexErrMsgTxt( "find. cluster  requires  one  or  two  input  arguments."); 

}  else  if  (nlhs  >  1)  { 

mexErrMsgTxt("find.cluster  has  one  output  argument."); 

} 

IN . rows  =  mxGetM (prhs [0] ) ; 

IN. cols  =  mxGetN (prhs [0] ) ; 

IN. ptr  =  mxGetPr(prhs [0]) ; 

if  (nrhs  ==  1)-C 
default. threshold  «  1; 
threshold  =  ftdefault. threshold; 

}  else  { 
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threshold  =  itixGetPr(prhs[l]) ; 

} 

OUT. rows  =  IN. rows; 

OUT. cols  =  IN. cols; 

plhsCO]  =  mxCr eat eFulK OUT. rows, OUT. cols,  REAL); 
OUT . ptr  =  mxGetPr (plhs [0] ) ; 

find„ clusters (OUT, IN, threshold) ; 
return ; 

} 
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